Nutrient-Loading Mitigation by Shellfish Aquaculture in Semi-Enclosed Estuaries

Bivalve aquaculture may provide a variety of ecosystem services including nitrogen extraction from estuaries, which are often subject to excess nutrient loading from various land activities, causing eutrophication. This nitrogen extraction may be affected by a combination of various non-linear interactions between the cultured organisms and the receiving ecosystem. The present study used a coupled hydro-biogeochemical model to examine the interactive effects of various factors on the degree of estuarine nutrient mitigation by farmed bivalves. These factors included bay geomorphology (leaky, restricted and choked systems), river size (small and large rivers leading to moderate (105.9 Mt N yr-1) and high (529.6 Mt N yr-1) nutrient discharges), bivalve species (blue mussel (Mytilus edulis) and eastern oyster (Crassostrea virginica)), farmed bivalve area (0, 10, 25 and 40% of estuarine surface area) and climate change (water temperature, sea level and precipitation reflecting either present or future (Horizon 2050) conditions). Model outputs indicated that bivalve culture was associated with the retention of nitrogen within estuaries, but that this alteration of nitrogen exchange between estuaries and the open ocean was not uniform across all tested variables and it depended on the nature of their interaction with the bivalves as well as their own dynamics. When nitrogen extraction resulting from harvest was factored in, however, bivalve culture was shown to provide a net nitrogen removal in the majority of the tested model scenarios. Mussels provided more nutrient mitigation than oysters, open systems were more resilient to change than closed ones, and mitigation potential was shown to generally increase with increasing bivalve biomass. Under projected future temperature conditions, nutrient mitigation from mussel farms was predicted to increase, while interactions with the oyster reproductive cycle led to both reduced harvested biomass and nutrient mitigation potential. This study presents the first quantification of the effects of various biological, physical, geomorphological and hydrodynamical processes on nutrient mitigation by bivalve aquaculture and will be critical in addressing questions related to eutrophication mitigation by bivalves and prediction of possible nutrient trading credits.


INTRODUCTION
Bivalve aquaculture is increasingly being acknowledged for the ecosystem services it provides beyond its primary goal of food provisioning (van der Schatte Olivier et al., 2018;Alleway et al., 2019;Smaal et al., 2019). Among these services, nitrogen extraction from near-shore ecosystems (e.g. estuaries, lagoons and bays) is of particular interest as these regions are often subject to excess nutrient loading from various land activities (e.g. agricultural runoff, waste-water treatment effluent) resulting in regional eutrophication (Cloern, 2001;Bricker et al., 2007). Owing to their capacity to clear particles from large volumes of water, bivalves have been proposed to mitigate these excess nutrients and be included in nitrogen-trading programs (Lindahl et al., 2005;Lindahl, 2011;Nielsen et al., 2016;Petersen et al., 2019). While bivalve harvest clearly contributes to a net extraction of nutrients, the net effect of wild or cultured bivalves on nutrient dynamics while present in the coastal environment is not as obvious. These organisms may potentially interact with several biogeochemical processes that regulate ecosystem functioning. Their direct consumption of food particles can lead to changes such as increased intra-and inter-specific competition for phytoplankton (Cloern, 2005), enhanced water clarity (Meeuwig et al., 1998;Newell, 2004), and alterations in plankton community size structure as a consequence of their size-selective retention efficiency (Strohmeier et al., 2012;Sonier et al., 2016). In addition, non-assimilated food can also affect ecosystem dynamics. Particulate wastes (e.g. faeces and pseudo-faeces) sink faster than original food particles (Callier et al., 2006), which may lead to a net transfer of material from the water column to the benthos. This process could result in particulate organic material retention within near-shore ecosystems (Cranford et al., 2007). Part of the assimilated food is also excreted in dissolved organic and inorganic forms that can be more readily exported to the open ocean. These same dissolved wastes (e.g. dissolved organic nitrogen and ammonia) may also stimulate local planktonic primary production (Prins et al., 1995) and consequently increase the residence time of nutrients in the system. The net role played by cultured bivalves in material exchange between near-shore systems and the coastal ocean is difficult to predict a priori as it involves several non-linear interactions between the bivalves and the supporting ecosystem and can be influenced by various environmental, physical and farming conditions. Climate change has the potential to increase nutrient loading, through increased precipitation (Nazari-Sharabian et al., 2018) and to intensify eutrophication, through increased water temperature and its effects on many biogeochemical processes in coastal waters (Rabalais et al., 2009). The ultimate effect climate change could have on the role of bivalves in nutrient cycling through coastal ecosystems remains, however, highly uncertain as competing influences are at play. For instance, some bivalve species may grow faster at higher temperatures (Guyondet et al., 2015;Steeves et al., 2018). Temperatureassociated changes in primary production (Brown et al., 2010) and planktonic communities (Allen and Wolfe, 2013;Mackas et al., 2013), however, may alter bivalve food availability.
Temporal shifts in the seasonal temperature cycle may also affect the interaction between bivalves and the pelagic ecosystem by changing the phenology of planktonic species and/or the phenology of bivalves, especially their reproductive cycle (Filgueira et al., 2015).
Intrinsic features of near-shore regions may also greatly influence their response to pressures such as nutrient loading, bivalve aquaculture and climate change. Specifically, inlet morphology plays a major part in controlling the intensity of exchange with the open ocean and, consequently, strongly influences the prevailing internal conditions and the sensitivity to external stressors (Filgueira et al., 2013;Panda et al., 2013;Roselli et al., 2013). The objective of the present study was to determine the overall net nitrogen retention versus extraction of bivalve aquaculture in semi-enclosed estuaries and explore some of the driving factors of this relationship. Specifically, we focused on the exchange between various types of embayments and the open ocean under various levels of nutrient loading, different shellfish culture species/biomass and climate-change-induced temperature increase as drivers of coastal ecosystem response. Numerical modelling provides an adequate framework for this investigation involving the non-linear intricacies of bivalveecosystem interactions, as suggested in previous research using coupled hydro-biogeochemical model simulations . This modelling could help to characterize the eutrophication mitigation potential of bivalve culture in the context of different environmental and biological scenarios relevant to coastal management.

Coupled Hydro-Biogeochemical Model
The RMA10-11 suite of models was used for this work (King, 1982;King, 2003) in a 2D depth-averaged configuration with a spatial resolution of 100 -250 m (distance between nodes). This same modelling framework was successfully applied to bivalveecosystem interaction studies in the past (Guyondet et al., 2010;Guyondet et al., 2015). Details of the model structure are given in these previous reports and the latest configuration, used in the present study, was similar to Filgueira et al. (2016) and provided in Appendix 1. Briefly, the model simulates the nitrogen cycle through a coastal pelagic ecosystem consisting of dissolved inorganic nitrogen (DIN), phytoplankton (Phyto), zooplankton (Zoo) and organic detritus (D). Cultured bivalves are fully integrated within this pelagic structure and their physiology including feeding, somatic growth and reproductive cycle is simulated using a dynamic energy budget (Kooijman, 2010) with parametrizations from (Filgueira et al., 2011;Filgueira et al., 2014) for mussels and oysters, respectively. The model also accounts for the effect of water temperature on all simulated biological processes. The benthic nitrogen cycle is not explicitly simulated, rather temperature-dependent rates of DIN efflux from the sediment are prescribed outside and inside commercial bivalve farms according to observations made in a coastal lagoon of the Magdalen Islands, Atlantic Canada (Richard et al., 2007). River, atmospheric, and open-boundary conditions were constructed from observations in the southern Gulf of Saint Lawrence (GSL), Atlantic Canada (Guyondet et al., 2015) and are characterized by weak mixed diurnal and semidiurnal tides (0.15 -0.60 m amplitude). Further forcing details are provided in Appendix 1.

Model Factors: Bay Geomorphology
In order to draw general conclusions, we chose to carry out numerical experiments on geomorphologically-idealized coastal systems, which are representative of shallow temperate estuaries (model mean depth of 4.6 m). Three geomorphologies were tested that only differed by their inlet configuration, leading to leaky, restricted and choked systems in increasing order of tide and exchange attenuation (Kjerfve, 1986). This strategy allowed spanning over a wide range of systems that real estuaries could be compared to and provided a full control on the degree of exchange between the near-shore area and the open ocean as reflected by the differences in water renewal time among the three systems (Figure 1), computed following Koutitonsky et al. (2004).

Model Factors: Environmental Drivers
A factorial design was developed to test the effects of several factors and their interactions on the response of the coastal ecosystem. In addition to the system's geomorphology, the following factors were tested. The first was nutrient loading, with two levels, where nutrient inputs were obtained from small and large river discharges resulting in moderate (105.9 Mt N yr -1 ) and high (529.6 Mt N yr -1 ) loading, respectively (Coffin et al., 2018). This factor also influences the overall water exchange between the embayment and the open ocean. The second factor was bivalve species, with two levels, being either the blue mussel (Mytilus edulis) or the eastern oyster (Crassostrea virginica). These species are two of the most commonly cultured bivalves in Atlantic Canada. The third factor, with four levels, was the surface area of the embayment covered by bivalve aquaculture (no aquaculture and low, medium and high aquaculture), which was based on the farm surface area coverage of the embayment (0, 10, 25 and 40%, respectively, Figure 1). The "no aquaculture" scenario (stocking coverage = 0%) served as a reference to test the effects induced by the presence of the bivalve farms. The density of bivalves in farmed areas was typical of culture practices in southern GSL with 94.1 mussels m -2 and 30.0 oysters m -2 (Drapeau et al., 2006;Comeau, 2013). The remaining three factors were related to climate change, with two levels each. Water temperature, sea level and precipitation reflected either present or future (Horizon 2050) conditions. Current conditions were obtained from observations in the southern GSL (Guyondet et al., 2015) while future conditions were derived from predictions at the regional scale based on observed trends (Vasseur and Catto, 2008;Taboada and Anadoń, 2012) and corresponding to an annual mean 2°C increase in temperature [with warmer summer peak and earlier Spring Appendix 1)], a FIGURE 1 | Inlet configuration and water renewal time for the three bay geomorphology types (left) and bivalve farm stocking areas for the three aquaculture scenarios (right). 0.5 m sea level rise and a 10% precipitation increase, forced through increased river discharge.
The combination of all the factors and their respective levels led to a total of 336 scenarios, which were simulated over a full year from June to June, preceded by a 90-day spin-up period in order for the system to reach an equilibrium state. For the purpose of the present study, exchange fluxes of all variables across an inlet section were recorded and then integrated over the year to estimate their net annual import/export balance. These net figures for all variables were also combined into a single estimate of the overall net import/export balance of the coastal ecosystems using the conversion coefficients reported in Table 1. Bivalve production was computed as the difference in individual biomass between the end and beginning of the simulation, scaled up to the total cultured population of each scenario. Among the climate factors tested, Filgueira et al. (2016) noted that temperature had by far the largest influence on the ecosystem response. Therefore, temperature is the only climate factor explicitly discussed in this study.

Eutrophication Mitigation Potential
As mentioned in the Introduction, some of the feedbacks from cultured bivalves to the ecosystem can lead to retention of inorganic/organic material within the near-shore area (i.e. estuary or bay). This retention was quantified by comparing the net annual exchange of material for a given aquaculture scenario with the exchange for the corresponding scenario without aquaculture. The eutrophication mitigation potential is defined here as the difference between the bivalve production, representing the material that can actually be extracted through harvest, and the retention of material defined above. Bivalve production had to be estimated at a fixed date (end of simulation), which makes the evaluation of the mitigation potential sensitivity to the various factors somewhat tributary to the phenology of the biological processes, in particular the bivalve reproductive cycle (see Discussion).

Statistics
Statistical analyses were conducted using R version 4.0.2 (R Core Team, 2020-06-22) operating in the RStudio version 1.2.1335 environment. Significance for all statistical tests was set at p < 0.05.
Our first goal with the application of statistics was to understand how three independent variables of interest (bivalve area coverage, bay morphology and species) interact to predict changes in response variables, namely primary production (PP) and export to the open sea of phytoplankton biomass (Phyto), detritus biomass (Detritus), dissolved inorganic nitrogen (DIN), zooplankton biomass (Zoo), and total nitrogen (TN). Response variables were expressed as a percentage change when compared to non-aquaculture scenarios. In order to meet normality assumptions prior to analyses, all response variables were subjected to Yeo-Johnson transformations (Yeo and Johnson, 2000) using the recipe() function in the recipes package (Kuhn and Wickham, 2022). QQ plots were produced to visually assess the normality assumption for all groups together and also for each group separately. Residuals versus fits plots were produced to check the homogeneity of variances. The absence of extreme outliers was confirmed using the identify_outliers() function from the rstatix package (Kassambara, 2021). Finally, three-way ANOVAs were performed using the Anova() function from the car package (Fox and Weisberg, 2019).
Our second statistical objective focused on the extractive capacity of bivalve aquaculture. Our null hypothesis was that the annual amount of nitrogen extracted from the system due to bivalve harvest was equal to the annual amount of nitrogen further retained within the system (export reduction) due to the presence of bivalve farms. A comparison of scenarios 'with' and 'without' cultured bivalves revealed that differences between paired scenarios were not distributed normally, nor symmetrically around the median. Hence the paired samples sign test was used to determine whether differences were consistent across scenarios. The same approach was utilized for assessing whether differences were consistent across species.
All statistical results, raw data, and annotated R-script can be found in the Supplementary Material (Appendices 2-4).

Embayments Without Aquaculture
As expected, all coastal system configurations tested in this study under "natural", i.e. without bivalve farms, and both present and future conditions, showed a net export of material (all forms of N combined) to the open sea ( Figure 2). The magnitude of this export is mainly driven by the inlet morphology and the river dischargethe more open the system and the larger the river, the greater the nitrogen exports to the open sea. The total N export ranged from 13.1 to 316.7 Mt N yr -1 for small rivers and 27.5 to 532.5 Mt N yr -1 for large rivers, while these systems received a total N river loading of 105.9 Mt N yr -1 and 529.6 Mt N yr -1 , respectively, of which 60% was of inorganic form (DIN) and 40% of organic form (Detritus) for both river sizes.
There are some obvious associations that can be made between nutrients and primary biomass/production within bays ( Table 2). High nutrient loading from rivers are associated with high concentration of all pelagic variables (Phyto, Zoo and DIN), and also with elevated rates of PP. With respect to bay morphology, there is a tendency for

Embayments With Aquaculture
For all scenarios, the addition of bivalve farms consistently led to a decrease in total N export to the open sea ( Figure 3A), signaling the presence of aquaculture-related retention mechanisms at the bay scale. Export reduction was influenced by a two-way interaction between stocking area and bay morphology, F(4, 270) = 6.53, p < 0.001, h 2 p = 0.09 (Appendix 2; the three way interaction including the species factor is non-significant). Stocking area has a highly significant effect on export reduction (p < 0.001) across all bay morphologies, with a strikingly large effect in choked systems (h 2 p = 0.54 for mussels, 0.53 for oysters). Bay morphology also has a highly significant (p < 0.001) and large (h 2 p = 0.13 to 0.17) effect on export reduction, although only at high and medium stocking area levels; at low stocking area levels, the effects of bay morphology on export reduction are either non-significant (mussels, p = 0.09) or significant but small in magnitude (oysters, p < 0.01, h 2 p = 0.04). The addition of bivalve farms increases grazing pressure (topdown control), which leads to a decrease in primary production in most scenarios ( Figure 3B). The low stocking area results show however that there is potential for a positive feedback from bivalves on primary production through the excretion of DIN (bottom-up control). As with total N export, PP is influenced by a two-way interaction between stocking area and bay morphology, F(4, 270) = 4.45, p < 0.01, h 2 p = 0.06; however, PP is also regulated by a two-way interaction between stocking area and species, F(4, 270) = 4.29, p < 0.05, h 2 p = 0.03. In more detail, the effects of stocking area on PP are significant (p < 0.05) for all bay morphologies, with the exception of mussels in leaky systems (p = 0.16). Where significant the magnitude of the stocking area effect ranged from a low effect size (h 2 p = 0.03 for oysters in leaky bays) to a large effect size (h 2 p = 0.17 for oysters in choked bays). Overall, these factorial ANOVAs revealed that both mussels and oysters can decrease PP and reduce total N exports to the open sea, particularly when cultivation activities are carried out at high intensity in choked systems.
The alteration of material exchange induced by bivalve farming is not uniform among pelagic components (Figure 4). Despite being the preferred bivalve food source, Phyto export was among the least affected, while Zoo showed the highest relative change. ANOVAs corroborated the data's inference that stocking area coverage increased retention across all pelagic components (Appendix 3). However, the nature of the stocking area effect varied across components. A two-way interaction between stocking area and bay morphology   While bivalve aquaculture reduced exports to the open sea, it may also provide an extraction service due to harvests. It was found that the magnitude of these contrasting effects differs between cultured species ( Figure 5). Export reduction mechanisms were stronger for oysters than for mussels (Appendix 4; p < 0.05, pairedsamples sign test), with the exception of leaky systems at low stocking density where there was no species effect (p = 0.077). By contrast, mussel production provides a much higher extractive potential than oyster production (p < 0.001, paired-samples sign test). Consequently, the mitigation potential was stronger for mussel farms (43.8 ± 5.6 -148.0 ± 35.1 Mt N yr -1 ) compared to oyster farms (1.1 ± 7.6 -29.4 ± 2.2 Mt N yr -1 ). The mitigation potential was highly significant across all culture scenarios (p < 0.001, pairedsamples sign test), with the exception of oysters cultured under high densities in leaky systems, for which there was no significant (positive or negative) mitigation effect. Net retention (negative mitigation) only occurred in a small subset of scenarios (15 out of 288), corresponding to oyster culture in leaky systems at medium or high stocking area coverage under future temperature conditions. Water temperature interacted with other factors such as bivalve species, stocking surface area and system morphology ( Figure 6). The mitigation potential of mussels increased with temperature for the majority of scenarios (54/72 points above the 1:1 line). In contrast, the mitigation by oysters tended to decrease at higher temperature and all cases of negative mitigation potential (i.e. net N retention) occurred in leaky systems and future temperature conditions. This result, although counterintuitive and also contradicting the faster growth of oysters at higher temperature (not shown, but see companion paper Filgueira et al., 2016), is caused by the strong influence of temperature on the reproductive cycle of oysters, which leads B A FIGURE 3 | Mean ( ± SD) relative changes in (A) total nitrogen export to the coastal ocean and (B) primary production, following the introduction of various cultured bivalve species at various stocking surface areas in the three bay types. Error bars represent variability introduced by other factors (i.e. nutrient loading and climate).
to an increased reproductive output (larger spawning), which in turn may decrease the oyster final biomass and production estimate. The positive relationship between mitigation potential and stocking area for a given morphology is preserved in future temperature scenarios except for oyster farms in leaky systems. The general rule stating that the stronger the flushing the more resilient the system is to changes (points closer to the 1:1 line) is verified for the temperature effect on mitigation potential of mussel farms, leading more enclosed systems to benefit more from mitigation in future climate scenarios. However, mitigation by oyster farms departs from this rule due to the influence of temperature on their reproductive cycle mentioned above. Finally, the clearer temperature response pattern provided by the mussel farm scenarios allows the viewing of a different behaviour between the small and large river systems (only visible for medium and high stocking areas). In systems with large rivers, higher flushing from wider openings tends to lead to a lower mitigation potential, consistent with the higher resilience hypothesis, while the reverse occurs in small river systems, especially at current temperature forcing. Moreover, while increased temperature tends to exacerbate the flushing effect in large river systems (stretch along the y-axis), the range of mitigation potential in small river systems becomes much narrower among flushing scenarios in future conditions (compressed y-axis range).

DISCUSSION
The use of idealized systems prevents the direct extrapolation of our results to existing estuaries, however the general conclusions drawn from these results and discussed in the following sections can be transposed to real systems with similar dimensions, loading, and exchange with the open ocean. The model outputs obtained in this study indicate that for all stocking area scenarios bivalve culture was associated with further retention of material within bays. This alteration of the bay exchange was not uniform across pelagic variables as it depended on the nature of their interaction with bivalves as well as their own dynamics. When the material extraction from harvest was factored in, however, bivalve culture was shown to provide a net nitrogen removal in the majority of scenarios tested. In addition, for the husbandry conditions simulated, mussel rather than  oyster farming provided the strongest potential for nutrient loading mitigation. Our results also support the idea that more open systems are more resilient to change. The corollary stating that more enclosed systems should benefit more from bivalve culture for nutrient mitigation, only holds for oyster farming scenarios. Moreover, the mitigation potential was shown to generally increase with bivalve stocking area. In future temperature conditions, the mitigation from mussel farms was predicted to increase, while interactions with the oyster reproductive cycle led to reduced harvested biomass and nutrient mitigation potential.

Simulated Embayments
The idealized embayments tested are representative of shallow estuarine systems in temperate regions with a potential for moderately to highly eutrophic conditions, with nutrient loading ranging from 21 to 106 kgN estuary ha -1 yr -1 , depending on the river size (Coffin et al., 2018). As can be expected for such river-driven systems, as opposed to those driven by coastal upwelling, the balance of material exchange at the ocean boundary results in a net export. The more open the system is to the ocean, the higher its export capacity becomes. In the moderate loading (small river) scenarios choked systems could only export 1/10 of the river input, leaky systems being able to export an amount equivalent to three times that same river loading. The apparent imbalance between inputs and export of these leaky scenarios is explained by the ability of the more open systems to also export part of the nitrogen provided by benthic loading. However, the exporting capacity seems limited as shown by the higher river loading scenarios where choked systems could only export 1/20 of the nitrogen received and leaky systems exported the equivalent of only one time what they received from river input. Hence, material retention within the estuary increases in more enclosed systems and for larger river loading. This material retention has been referred to as the filtering capacity of an estuary and for nitrogen has been reported to range from 0 to over 50% of the river loading on an annual basis (Arndt et al., 2009 and references therein). Benthic processes such as burial and denitrification may contribute to the net removal of this retained material (Seitzinger et al., 2006), however, they were not explicitly detailed in the model structure so their contribution and change among model scenarios was not evaluated. In real ecosystems, the retained inorganic/organic material can fuel the internal productivity of higher trophic levels, both pelagic and benthic, as exemplified by the higher Zoo predicted in this study for enclosed embayments. Counter-intuitively, these systems show a smaller PP rate than more open ones ( Table 2). This result must be regarded from an efficiency perspective as it relates to the importance of the balance between the biological and hydrodynamic time-scales of a coastal system (Officer, 1980;   Kimmerer et al., 1993) that dictates how much of the nutrient loading can be used by local biological processes . In more open bays, the faster exchange dynamics are more favorable to the fast processes of phytoplankton production (Eyre and Twigg, 1997), while in enclosed or restrictedopening bays, slower exchange leads to a more efficient transfer of energy towards secondary producers (i.e. zooplankton). Part of the reduced pelagic PP in the latter can also be explained by the overall higher concentration of suspended particles, which can increase light attenuation and potentially decrease PP rates, even in bays with ample nutrient supply (Meeuwig et al., 1998). Moreover, in shallow estuaries, high nutrient loading leading to eutrophic conditions often manifests itself through the proliferation of opportunistic macroalgae, e.g. Ulva spp. (Valiela et al., 1997). Accounting for this other type of primary producer would provide a more realistic representation of these system dynamics (Lavaud et al., 2020) as they may affect the filtering capacity of the bay by being likely less easily exported than phytoplankton and by storing large nitrogen pools that would not be directly amenable to mitigation through bivalve farming. Climate change, leading to warmer future temperatures, could possibly increase metabolism in coastal systems (Angilletta et al., 2002;Brown et al., 2004), resulting in faster biological dynamics, without directly affecting the hydrodynamic time-scale. Consequently, these future conditions could favour an increase in the filtering capacity or material retention of these nearshore regions. However, warming could also alter coastal ecosystems more profoundly as it might not affect all communities and trophic levels in the same way or intensity (Carr and Bruno, 2013;Mertens et al., 2015;Alexander et al., 2016;Ullah, 2018). Hence, future changes in the filtering capacity of estuaries and coastal systems remain highly uncertain and warrant further research at the ecosystem scale.

Bivalve Influences
Cultured bivalves feed preferentially on phytoplankton and consequently build their tightest link with this component of the ecosystem (Filgueira et al., 2011;Saraiva et al., 2012). The net balance of this interaction results, however, from the competition between the top-down feeding pressure and the bottom-up control exerted by the ammonia excretion of these mollusks (Prins et al., 1995;Smaal et al., 2013). The positive response of PP in some of the low stocking area scenarios is a direct manifestation of this bottom-up feedback. Nutrient loading also contributes to the bottom-up side of this balance, with higher loading generally counter-acting more of the feeding pressure from cultured bivalves. Finally, this balance is modulated by the system's morphology and temperature as can be partially seen in Figure 3B (see Filgueira et al. (2016) for a detailed analysis).
When considering the impact of bivalve aquaculture on the exchange of material with the open sea, phytoplankton appears as the most resilient variable. The low impact on it likely results from its fast turn-over rates and the feedback mechanism mentioned above. Conversely, zooplankton exchange is the most affected due to its slower dynamics, the competition for food with the bivalves (Gianasi et al in review; Granados et al., 2017a;Granados et al., 2017b) and, although marginal, the direct bivalve predation on zooplankton (Pace et al., 1998;Davenport et al., 2000;Trottet et al., 2008;Peharda et al., 2012). Despite the overall low contribution of zooplankton, representing only 1.8 to 3.9% of total N exported in all scenarios tested, this last result indicates that the bivalve aquaculture signal can be amplified through the trophic levels both in and out of these nearshore systems. Similarly, the retention efficiency of bivalves as a function of food particle size (Riisgard, 1988;Strohmeier et al., 2012;Sonier et al., 2016), coupled with the strong feeding pressure exerted in aquaculture settings, could affect the pelagic community size distribution at a system scale. According to size-spectrum theory (Sheldon et al., 1972;Andersen et al., 2016), bivalve culture could then affect the way energy propagates through trophic levels. Further work is required to better understand how and how much these cascading effects might change coastal ecosystem functioning and the services they provide, and to examine the balance between these negative effects and the positive effects that accrue from added food production and nutrient loading mitigation from bivalve aquaculture.
All tested scenarios showed an increase in estuarine filtering capacity (i.e. material retention) in the presence of bivalve culture farms. This increment in retention ranged from 0.3 to 4.3% and from 1.5 to 18.5% of river N loading for large and small rivers, respectively. Despite the lower stocking density typically used in oyster versus mussel farming (30 vs 94 ind. m -2 of farm in the present model applications) and the slightly lower individual feeding rate resulting from the mean size and environmental conditions over the simulated period, oyster culture provided a slightly higher increase in retention compared to mussel culture for the same farm coverage. This species difference rests mostly on a feeding efficiency variance at the individual level and in particular a higher feedback by mussels through excretion, which gets amplified by the difference in stocking density. Finally, in absolute value, material retention from bivalve farms is larger in systems with stronger flushing intensity, which can be explained by the retention/export being already high/low in more enclosed systems, thus limiting the potential for further retention from bivalve farms. This result is also in accordance with the negative relation between the fraction of loading exported and freshwater residence time reported by Dettmann (2001).

Mitigation Potential
When the material retention effect is examined in the context of the net material removal provided by the harvest of bivalve production, a positive nutrient loading mitigation potential is predicted for all mussel and most oyster culture scenarios tested in this study. Moreover, such mitigation potential can reach important fractions of the river loading, with values ranging in high loading conditions (large river) from 8.86 to 36.54% and from 0.21 to 6.31% for mussels and oysters, respectively, and in moderate loading conditions (small river) from 35.14 to 132.69% and from 1.57 to 27.80%, respectively. The stronger mitigation potential provided by mussel farming can be attributed, at the individual level, to the faster growth of this species (12 -18 months and four years to reach harvest size for mussels and oysters, respectively) and, at bay-scale, to the higher densities they are typically reared in (94 mussels vs 30 oysters per square meter of farm in the present study). The higher mitigation efficiency predicted in lower loading conditions reflects the limits to the assimilative capacity that bivalve farms can provide.
The tested mussel culture scenarios show an interesting interaction between loading and flushing intensities. The reversal of the relationship between mitigation potential and flushing, depending on loading intensity, highlights the importance of the absolute loading level for the vulnerability of coastal ecosystems to eutrophication and for the mitigation benefit cultured bivalves can provide, as well as the tight link between hydro-and biological dynamics in these regions (Officer, 1980;Kimmerer et al., 1993;Dame et al., 2000). From this result, it appears that a loading inflection point exists beyond which more flushing becomes detrimental to the mitigation potential of bivalve farms. Further testing with incremental loading levels would be required to refine this relationship. In addition, spatial heterogeneity was not assessed in this study, given the idealized morphology of the test bays used. In natural systems, however, the realized mitigation could depend on the distribution of wild and cultured bivalves in relation to the nutrient loading outfall location and the distribution of water residence time (Gray et al., 2021). A spatially-explicit approach would then be critical in any real-life mitigation potential assessment.
Future warmer temperatures could have competing influences on mitigation potential. As stated earlier, a general increase in biological activity is expected that could lead to increased material retention. At the same time, cultured bivalves could also benefit from increased metabolism to reach higher production rates, which would favour nutrient removal. However, the actual future response of bivalve production in a given region will depend on how future temperatures relate to each species' tolerance range (Steeves et al., 2018), how tightly adapted they are to current conditions and how quickly they can evolve in response to these changes (Pernet et al., 2007;Thomas and Bacher, 2018;Thyrring et al., 2019). According to the scenarios tested, not accounting for potential adaptation effects, mussel culture could provide increased mitigation at these future temperatures. The picture provided by the oyster culture scenarios is much less clear and even tends to indicate movement in the opposite direction despite the tolerance for higher temperatures of that species. Part of this somewhat surprising result can be explained by a method artefact, which comes from the necessity in our modelling framework to evaluate bivalve production at a fixed date at the end of the simulation period and how this can interact with the changes in the phenology of the reproductive cycle of the species considered. The potential for oyster production is actually higher at higher temperature, but the realized production is sometimes lower because of the evaluation at a fixed date. However, this result highlights the importance of the reproductive cycle in optimizing bivalve culture for nutrient loading mitigation as these species can invest large fractions of their total tissue weight in gametes [over 40% for American oysters (Choi et al., 1993) and similar values (35 -38%) for blue mussels (Sukhotin and Flyachinskaya, 2009;Hennebicq et al., 2013). Such phenological change considerations tie to the broader seasonality of the various processes involved in the filtering capacity of nearshore systems (Brion et al., 2008) and ultimately in the nutrient loading mitigation potential from bivalve culture, especially in temperate regions.
Bivalve culture has been proposed as a mitigation tool in nutrient trading schemes in Europe and the United States with studies evaluating its costs/benefits against more traditional methods (Lindahl et al., 2005;Lindahl, 2011), assessing the removal potential at farm-scale (Clements and Comeau, 2019;Bricker et al., 2020) and operationalizing and optimizing the production for this specific mitigation service (Petersen et al., 2014;Taylor et al., 2019). The results of the present study concur on the ability of bivalve farming to provide such a nitrogen extraction service and detail the influence of various drivers on the realization of this mitigation. In a management context, the seasonal considerations mentioned earlier and the tight and complex links between bivalves and their environment advocate for the use of a dynamic and spatially-explicit approach, such as the modelling tool used in this study, to provide an integrated assessment of the mitigation potential of bivalve culture within the broader coastal ecosystem functioning.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.  generated the data required to build the forcing conditions used in this modelling study.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2022. 909926/full#supplementary-material Seven supplementary files are provided: 1. Appendix 1: Model structure and forcing (equations, parameters and river, atmospheric and open ocean physical and biochemical forcing conditions). 2. Appendix 2: Commented R script for statistical analyses related to Figure 3 and corresponding data (2 files). 3. Appendix 3: Commented R script for statistical analyses related to Figure 4 and corresponding data (2 files). 4. Appendix 4: Commented R script for statistical analyses related to Figure 5 and corresponding data (2 files).