Production Characteristics and Optimization of Mitigation Mussel Culture

Taylor, Daniel; Saurel, Camille; Nielsen, Pernille; Petersen, Jens Kjerulf

doi:10.3389/fmars.2019.00698

ORIGINAL RESEARCH article

Front. Mar. Sci., 12 November 2019

Sec. Marine Fisheries, Aquaculture and Living Resources

Volume 6 - 2019 | https://doi.org/10.3389/fmars.2019.00698

Production Characteristics and Optimization of Mitigation Mussel Culture

$\r\nDaniel Taylor*$ Daniel Taylor^*

Camille Saurel

Pernille Nielsen $Jens Kjerulf Petersen\r\n$ Jens Kjerulf Petersen

Danish Shellfish Center, National Institute of Aquatic Resources, Technical University of Denmark, Nykøbing Mors, Denmark

Bivalve environmental services have become a focal point for their inherent role in the management of eutrophication, while active cultivation has become increasingly acknowledged as a mechanism for integrated nutrient reduction. In recent years, cultivation practices designed specifically for nutrient extraction have emerged; “mitigation culture.” While modeling efforts have been able to describe expanded potential of these services, only a single commercial pilot scale, real-world demonstration, has been documented. Over two production seasons (2017–2018), the optimization of nutrient extractive potential of mussels (Mytilus edulis) at full commercial-scale was evaluated by first testing multiple density configurations of conventional longline-spat collector setups and potential harvest times, then by comparing different cultivation technologies at three farms. Potential biomass volumes of 770–1700 t with longlines and 2100–2600 t on nets was demonstrated in full-scale production (18.8 ha), yielding 0.6–1.27 t N ha^–1 and 0.04–0.1 t P ha^–1, and 1.63–2.0 t N ha^–1 and 0.1–0.12 t P ha^–1 respectively. In general, 1 t of harvested mitigation mussels will yield 13.7 kg N and 0.9 kg P. Winter harvests exhibited higher yields (103–124%) than early spring harvests on optimized configurations, favoring an abbreviated production season. Production potential was similar between sites, despite differing environmental conditions, indicating eutrophic waters are suitable for expanded mitigation production. This study presents for the first-time production data of mitigation mussels utilizing different configurations and technologies to maximize yield and nutrient extraction potential.

Introduction

Global anthropogenic nutrient enrichment of coastal marine waters and estuaries is recognized as one of the principle negative drivers of environmental and ecological change in the biosphere (Rabalais et al., 2009; Doney, 2010; Jessen et al., 2015). The EU, and especially northern European waters, have experienced persistent and widespread marine eutrophication over the past century (Billen et al., 2011), which has given rise to legislation such as the Water Framework Directive (WFD) requiring European coastal waters to reach “Good Environmental Status” by 2020 (Borja et al., 2013). As of 2018, only 46% of all EU coastal waters had achieved this level or higher (Kristensen et al., 2018). Similar comprehensive legislation have been adopted in other industrialized parts of the world with persistent coastal eutrophication, such as the Clean Water Act in the United States (Copeland, 2012).

Due to primarily implemented modifications of wastewater treatment plants, point sources of nutrients are becoming less influential in the management of eutrophication, abatement efforts are focused on diffuse agricultural and atmospheric sources (Andersen et al., 2014). In Denmark > 60% of the land area is currently used for agriculture. General regulatory regimes have effectively reduced loads into the aquatic environment, however, further intensification of these reductions (mitigation) are anticipated to be expensive (Hasler et al., 2015). Spatially-explicit targeted modifications to regulatory mechanisms and related abatement mechanisms is required to further progress toward WFD goals for coastal waters, and is expected to be a challenge (Dalgaard et al., 2014; Konrad et al., 2014; Hashemi et al., 2018).

Variegated efficacy in nutrient management programs, persistent internal nutrient loads, and atmospheric deposition has elicited opportunities for the management of nutrients within the marine environment. In this context, much of the literature has focused on lower trophic assemblages.

Mussel cultivation has been proposed as a mitigation measure (Haamer, 1996; Lindahl et al., 2005; Gren et al., 2009; Stybel et al., 2009; Schernewski et al., 2012; Petersen et al., 2014, 2016, 2019a; Nielsen et al., 2016; Galimany et al., 2017). So-called “Mitigation Mussel Cultivation” (with Mytilus edulis) has previously been described in conceptual terms and only field tested at one commercial scale farm in Denmark over a single growing season (Petersen et al., 2014; Nielsen et al., 2016).

The mitigation concept is framed by relatively simple mass balance principles. Inorganic nutrients introduced into the marine environment are assimilated by phytoplankton and bacterioplankton. A large proportion of this planktonic mass is immobilized when filtered by suspended mussels. When mussels are harvested, assimilated nutrients are removed from the local marine environment. This mass balance principle thus focuses on the processes and subsequent nutrient extraction potential that are captured through harvest. In contrast to cultivation modes optimized for the provision of mussels for human consumption, mitigation culture aims to maximize the nutrient content of a given cultivation area. While cultivation for human consumption also generally intends to maximize biomass yield per farm, size, condition, and quality of the mussels are equally as important as total yield (Pérez-Camacho et al., 2013); so cultivation durations are longer and may incorporate intermediate steps to facilitate uniformity in size (i.e., thinning, socking). While production characteristics have been described for cultivation practices designed for human consumption, this has not yet been documented for mitigation production, where mussel size can be heterogeneous.

In the present study, differences in production yields of mitigation mussels were examined at three different farms in conjunction with modified conventional practices and alternative technologies for an optimized mitigation production. To further explore nutrient potential, a farm-scale model was used to identify potential production limitations. We hypothesize that increasing density and depth of conventional spat collectors, and the use of technologies with greater total surface area would proportionally increase nutrient extractive potential in eutrophic waters.

Materials and Methods

Study Areas

Limfjorden is a shallow (mean depth = 4.8 m), microtidal, partially stratified system of inland sounds with an approximate mean volume of 7100 10⁶ m³ (Hofmeister et al., 2009; Figure 1). A salinity gradient is formed along the west-east axis as the salinity in the Kattegat is diluted by freshwater inputs from the Baltic Sea (∼22–33 PSU) compared the high-saline water inflow from the North Sea (>32 PSU). The larger Limfjorden catchment area is composed of 90 subcatchments, a total area of 7528 km² of which ∼62% is agricultural. Drainage inputs approximate 2700 10⁶ m³ per year, with estimated annual loads to the estuary of 11,000–17,800 t N and 220–400 t P (Thodsen et al., 2018). Large areas of the Limfjorden experiences episodic hypoxia over the year (Møhlenberg, 1999; Conley et al., 2011), and the whole fjord system is, according to the WFD, classified as in “poor ecological condition” (Miljø-og Fødevareministeriet, 2016).

FIGURE 1

Figure 1. Map of mainland Denmark and location of three study sites, inset to a map of Europe. Limfjorden catchment is illustrated in dark gray. SKIV, Skive Fjord site; DV, Dråby Vig site; SALL, Sallingsund site.

Three sites were employed in this study (Figure 1): Skive Fjord (SKIV), Dråby Vig (DV), and Sallingsund (SALL), each with an area of 21.87, 23.06, and 14.76 ha respectively. These sites are situated in three hydrographically distinct areas within the greater estuarine system. SKIV is located centrally within Skive Fjord and at the junction of Lovns Bredning (Lovns Broad) and receives freshwater from Karup Å (Karup River) and Hjarbaek Fjord. Skive Fjord is the most eutrophic part of Limfjorden and can be characterized as a separate water body with its own (sub)catchment (Carstensen et al., 2013), with annual chlorophyll-a concentrations averaging 16.1 μg l^–1 and up to 58 μg l^–1 (NOVANA ODA 2018)¹. DV is positioned in Løgstør Bredning (Løgstør Broad) along the western extent of a large basin in mid-Limfjorden. SALL is located at the northern extent of Sallingsund (Salling Sound), is relatively narrow and experiences higher bidirectional current velocities than the other two sites. Water depths in the system are predominantly influenced by wind driven influx of seawater from the North Sea and water height differences between the North Sea and Kattegat. Site depths range between 5 and 7 m, based on bathymetry across the farm area and water height.

Cultivation and Test Configurations

Three preexisting commercial mussel farms at the study sites were used to test configurations and substrates. In 2017, standard longlines were used to test spat collector configurations. Two full-scale mitigation farms were utilized in May 2017, with each 18 standard long lines for experimentation and 72 for normal production (non-experimental): one farm in SKIV, one in DV. Three test lines were deployed within a commercial-scale multi-use mussel farm in SALL in May 2017. All standard long lines were 148 m in length, equipped with 5 cm polypropylene belt spat collectors, buoyed with 20 L plastic buoys, and counterweighted with concrete blocks at an interval of ∼5 m. Two factors were tested at each site to evaluate linear yield in 2017: spat collector spacing (40 cm separation between loops, with 30 or 60 cm loop width as tied on the mainline) and spat collector length (2 or 3 m depth); treatments are labeled 30 2, 30 3, 60 2, and 60 3 m (Figure 2A). Industry standard for spat collection for food production is 60 2 m. Treatments were crossed and placed by stratified pseudorandom assignment to account for spatial effects (e.g., northern half vs. southern half). Ice coverage in early 2018 (20/2–20/3) necessitated reduction in buoyancy and lowering of lines below ice formation (1–2 m) at all sites. Lines were raised again after ice retreat with reapplication of buoyancy and harvested soon after.

FIGURE 2

Figure 2. (A) Schema of the belt spat collector treatments and example of the substrate material (left). Loop spacing on the mainline was tested at 30 and 60 cm between loops; loop length from the mainline was tested at 2 and 3 m. Loop width is set at 40 cm. Treatments were established to fill a complete line (SKIV, SALL) or split on a line so that one line had two treatments (DV). Spat collectors are suspended from a mainline with 20 L buoys and concrete blocks (not illustrated here) for buoyancy maintenance. The water surface and sediment surface are represented in a 5–6 m deep water column. (B) Schematic of a 17.5 cm × 17.5 cm tube-net. The net is suspended from an air-filled polyethylene tube, and a head rope.

In 2018, at the same sites, a series of tube-nets were tested in conjunction with spat collectors at commercial scale within the same sections of each farm used in the 2017 cycle, i.e., the full farms were exploited. In SKIV, twenty 100 × 3 m Smartfarm² tube-net systems (Figure 2B) of alternating mesh size (ten tube-nets of 17.5 × 17.5 cm, ten of 25 × 28 cm) were deployed along with four longlines equipped with belt spat collectors at 30 cm spacing and 2 m depth (30 2 m). In DV, four 100 × 3 m tube-nets of alternating mesh size (17.5 × 17.5 cm, 25 × 28 cm) and four longlines with belt spat collectors at 30 cm spacing and 2 m depth (30 2 m) were deployed. One 100 m tube with three 30 × 3 m nets of differing mesh size (17.5 × 17.5 cm, 25 × 28 cm, and 20 × 20 cm) was deployed in SALL; along with two parallel longlines coupled by 200 L Xplora buoys each equipped with rope ladders (30 × 15 cm grid size) at 1.5 m depth and 1 m attachment spacing. Alternative treatments are labeled accordingly: Net 17.5, Net 20, Net 25, and Ladders. Additionally, one full longline was deployed with belt spat collectors at 30 cm spacing and 2 m depth (30 2 m) in SALL. Longline treatments in 2018 were selected based on experiences from the 2017 growth season.

Environmental Monitoring

All sites were sampled at approximately monthly intervals from June 2017 to the end of the year, then with less frequency until June 2018 when regular intervals were followed until the end of the year. At each event, a CTD cast (EIVA Arop) was performed at 0.1 m depth intervals to profile the water column for conductivity (salinity in PSU), temperature (°C), and chlorophyll-a (chl-a) fluorescence (μg l^–1). Secchi depth was recorded from the non-shaded side of the boat in triplicate in the approximate center of the farm (in between lines) and directly outside of the farm (approximately 300 m from the farm, at the continuous monitoring station detailed below). Water samples were captured at 2 m depth with a 3 L Ruttner sampler in triplicate for determination of chl-a. Chlorophyl-a samples were filtered through 25 mm Wattman GF/F filters, filters transferred to 10 ml of 96% ethanol, covered and light-sealed in refrigeration overnight, then measured for fluorescence before and after acidification in a Turner 10AU fluorometer (Holm-Hansen et al., 1965).

Continuous monitoring stations were established and moored approximately 300 m outside of each farm. Temperature and salinity were monitored continuously at all sites throughout the growing seasons (June 2017–December 2018). Chlorophyll-a and turbidity were monitored continuously from July 2017 to December of 2017, where after instrumentation was removed in anticipation of ice formation. Chlorophyll-a and turbidity monitoring was reestablished at all sites in June 2018 and remained until harvest in the winter. Buoys integrated with sensor capabilities were constructed and deployed at SKIV and DV. These buoys were equipped with a Star ODDI DST CTD at 5 min logging intervals, positioned at 2 m from the surface to track temperature (°C) and salinity (PSU). A collocated Turner Cyclops 7 chl-a fluorometric sensor and turbidity sensor were used in conjunction with a custom-built Arduino-based data logging system; 1 s readings were averaged into 30 s bins and recorded on an SD card. Sensors were serviced semi-monthly and calibrations performed with extracted samples as described above. In 2017 and until August of 2018, a Turner SCUFA was used to continuously (5 min frequency) record chl-a values at SALL, thereafter it was replaced with a Cyclops 7 logging buoy.

Biomass Sampling

Following mussel settlement in the late spring and early summer, regular sampling events were carried out to track mussel biomass and densities. Monthly bulk samples were captured in triplicate directly from belt spat collectors; a 30 cm section was removed from mid-depth of the spat collector loop (centered on 1 m from the upper extent of the loop). Single grid cells were excised from nets at upper- (0–1 m), mid- (1–2 m), and lower-depth (2–3 m) sections of the net by SCUBA diver. Samples were captured at regular spatial intervals in order to account for localized effects; e.g., at 1/3 along the line or net. Sampling locations were replicated for each treatment across each farm according to the randomized assignment and location within the farm (e.g., triplicate bulk samples from location X₁Y₁, X₁Y₂, X₂Y₁, X₂Y₂); sample replication is represented (as DF) in the results. Sampling repetition was elevated at harvest times to enhance estimations of biomass distribution within the farm. From these samples, total wet weight of the mussel biomass was determined after separation from byssus, detritus, and other intraspecifics. Subsampling of each replicate was performed to derive mean shell lengths (n = 100) and morphometrics (n = 10); mussels were selected haphazardly from the bulk sample.

Substrate length and surface area of each substrate type was calculated to evaluate yields in terms of substrate density. The ladders and nets are constructed of rope, which is approximately 1.4 cm in diameter, and have a simple cylindrical surface area of 476 cm² m^–1; while the belt spat collector exceeds 1000 cm² m^–1. These values were scaled to the line or net basis to standardize surface area biomass yield, based on total length of substrate per line or net.

All subsamples were kept frozen in sealed containers at −20°C until processing, which occurred no longer than 3 weeks following acquisition. Mussels were individually weighed whole, measured for shell length, width, and height; dissected and weighed for tissue and shell wet weight; dried at 80°C for > 3 days, weighed for dry tissue and shell weight. All individual weights were recorded to the nearest 0.1 mg. Condition Index (CI) factors are reported by the widely used allometric index: DW/SL³, where DW = dry tissue weight (mg) and SL = shell length (cm).

On harvest dates, mussels were stripped from spat collectors and conveyed into big-bags (∼2 m³), then weighed. A selected number of lines were harvested from SKIV (2 big-bags of each 18 lines, ∼20% of each line) and DV (6 total lines) in December 2017 to account for treatment replication; SALL was not harvested in the winter. Remaining 2017 lines in SKIV, DV, and SALL were harvested in late March and early April of 2018. In the 2018 growth season, tube-nets and spat collectors in SKIV were harvested in late November of 2018, while tube-nets were harvested at DV in January 2019. Samples were extracted during harvest from each treatment in an identical method to the regular biomass sampling program to verify estimation predictability of the method. A subset of harvest weights was additionally captured by measuring total harvested material per treatment by individual line. These weights were then regressed against sample estimates to assess sampling accuracy.

Subsamples of mussels from all treatment groups were haphazardly selected, dissected, and pooled for nutrient analyses by each farm and harvest date. Tissue samples were kept frozen in sealed containers at −20°C until processing. Samples were homogenized with a Krups Speedy Pro homogenizer and subsequently analyzed in triplicate per time period and farm for nitrogen by the Kjeldahl method (ISO, 2009) and total phosphorus by the spectrophotometric method (ISO, 1998).

Total potential harvest yields were calculated by scaling biomass per meter or square meter (dependent on substrate) to total substrate length within a “model farm.” A model farm is formed by 3 sets of 30 longlines or ladders (200 m length), or 80 tube-nets (120 m length), in 18.75 ha. Longline belt spat collectors have 1543, 1080, 1971, and 1380 m of substrate length for the 30 2, 60 2, 30 3, and 60 3 m respectively. Alternative technologies have 2160, 5143, 7200, and 8229 m length of substrate per line for ladders, Net 25, Net 20, and Net 17.5 respectively. Due to harvest equipment space requirements, tube-net systems are spaced at slightly greater intervals. Mass proportions of tissue, shells, and byssus captured in biomass and morphometric sampling were used to determine potential nitrogen and phosphorus extractive potential. Dry weight shell and byssus nutrient values were derived from previously reported values (Petersen et al., 2014): 0.98 g N and 0.0048 g P 100 g^–1 shell, and 11.02 g N and 0.077 g P 100 g^–1 byssus.

Statistical Analysis

Biomass was quantified on a basis of per-meter of substrate length for traditional spat collectors, and by square meter for nets. To compare yields by technology, as spat collectors were sampled in one dimension, while nets and ladders were sampled as grid cells, biomass yield was further transformed to substrate length; termed linear substrate yield. Biomass yield and morphometrics in conventional spat collector treatments (spacing and length) were analyzed by two-way ANOVA (one-way ANOVA in SALL due to loss of one treatment group), and significant difference between means was determined at a p ≤ 0.05. Alternative technologies were analyzed by one-way ANOVA, and significant difference between treatment means was determined at a p ≤ 0.05. Where significant differences were found, means were then compared by Tukey-Kramer’s Honest Significant Difference test. Analysis of areal yields (biomass t ha^–1) across farms was analyzed as a randomized complete block model with a standard least squares regression for all treatments (except Ladders and Net 20 as these were present at only one site) and harvest periods, while farm sites are random blocks. Data normality was assessed by the Shapiro–Wilk test, and equality of variance by the Bartlett test, where p > 0.05. Temperature, salinity, and chl-a were averaged by day, and compared between sites by the Kruskal-Wallis test (Helsel and Hirsch, 2002).

As elements within a farm are not spatially independent, analysis of biomass and morphometric factors in SKIV (2017 and 2018) and DV (2017) were additionally analyzed for spatially dependent variation to assess potential food depletion or other spatial phenomena influencing mussel growth. DV was not spatially analyzed in 2018 as the experimental section of the farm included only four nets and two long lines, which was below the power threshold for detecting spatial effects. Sampling positions were transformed to grid positions relative to line and farm dimensions. Directionality of yield and condition index was analyzed by linear regression. Central tendency of production and biometric data was assessed by quadratic polynomial regression (curve formation toward the center) and normalization of Euclidean distances (i.e., center of the farm = 1, edge = 0), thereafter testing trends by linear regression. Global spatial autocorrelation was tested by constructing inverse spatial weight matrices and calculating Moran’s I (Moran, 1950), a commonly-used inferential statistic in spatial ecology. For 2017 farms, where traditional spat collectors were used Moran’s I is determined by a 2-dimensional array, defined in Equation (1); while nets used a 3-dimensional array as samples were extracted at depth intervals in addition to latitudinal and longitudinal intervals, defined in Equation 2 (Marwan et al., 2007). The locations of clustered magnitudes were assessed by calculating Getis-Ord G_i^∗ as in Equation 3 (Ord and Getis, 1995), a commonly employed Local Indicator of Spatial Association in spatial statistics. Statistical analyses were conducted using JMP Version 13, except for spatial statistics, conducted in R (2017) and with the spdep package (Bivand et al., 2013).

I = \frac{N}{S_{0}} \frac{\sum_{j = 1}^{d_{1}^{*} d_{2}} \sum_{i = 1}^{d_{1}^{*} d_{2}} δ_{i j} (x_{j} - \bar{x}) (x_{i} - \bar{x})}{\sum_{i = 1}^{d_{1}^{*} d_{2}} {(x_{i} - \bar{x})}^{2}} (1)

I = \frac{N}{S_{0}} \frac{\sum_{j = 1}^{d_{1}^{*} d_{2}^{*} d_{3}} \sum_{i = 1}^{d_{1}^{*} d_{2}^{*} d_{3}} δ_{i j} (x_{j} - \bar{x}) (x_{i} - \bar{x})}{\sum_{i = 1}^{d_{1}^{*} d_{2}^{*} d_{3}} {(x_{i} - \bar{x})}^{2}} (2)

G_{i}^{*} = \frac{\sum_{j = 1}^{n} δ_{i j} x_{j} - \bar{x} \sum_{j = 1}^{n} δ_{i j}}{S \sqrt{\frac{[n \sum_{j = 1}^{n} δ_{i j}^{2} - {(\sum_{j = 1}^{n} δ_{i j})}^{2}]}{n - 1}}} (3)

Results

Environmental Monitoring

Time series of monitored parameters for each growing season are presented in Figure 3. As the sites are situated in the same large-scale watershed (Limfjorden) and are in relative close proximity to each other, temperature differences between sites were minimal. SALL and DV had higher salt concentrations (29.8 ± 1.9 and 28.0 ± 1.15 PSU respectively) than SKIV (24.03 ± 1.37 PSU) (X²₂, ₄₉₅ = 1058.11; p < 0.01). Chlorophyll-a concentrations varied between sites and within sites between years (X²₂, ₄₉₅ = 38.23; p < 0.01); SKIV generally exhibited higher concentrations over the growing season (June-October) than other sites (SKIV = 5.49 ± 2.76, DV = 3.56 ± 1.74, SALL = 3.58 ± 2.33 μg l^–1). From the discrete monitoring program, Secchi depth was on average 0.8 m deeper within the farms (3.8 m) compared to outside of the farms (3 m) and tended to be similar between all farms (ANOVA, p > 0.1).

FIGURE 3

Figure 3. Time series of environmental conditions at each site from continuous monitoring at 2 m depth. Left panel: orange = temperature (°C), blue = salinity (PSU). Right panel: green = chl-a (μg l^–1).

Biomass Growth and Yields

Biomass accretion between sites was similar [F₍₂, ₂₄₈₎ = 1,06, p = 0.3466], while inter-annual [t₍₄₃₎ = 1.90, P = 0.0633] and seasonal variability [t₍₇₃,₉₉₎ = −3.34, p = 0.001] affected yields to a greater extent. In general, 2 m collectors provided higher yields per meter collector [F₍₃, ₂₀₉₎ = 1.06, p < 0.001], and the 30 cm spacing configuration increased areal efficiency of yields compared to a conventional longline setup [F₍₃, ₂₀₉₎ = 30.93, p < 0.001]. Mussel size and condition generally followed areal density of mussels (Supplementary Table A2), which is related to substrate density and settlement. Nets yielded less biomass per meter substrate than belt spat collectors (Supplementary Table A1), and smaller mussels (Supplementary Table A2), but given higher total surface area, total yields were higher [Table 1 and Supplementary Table A4; F₍₇, ₃₄₃₎ = 86.91, p < 0.001]. Seasonal growth of each treatment is represented by per-meter substrate biomass yield in time, by site, in Figures 4A,B. Summary biomass yields, and densities are reported for all sites and by treatment in Supplementary Table A1 and summary morphometrics in Supplementary Table A2 with statistics. Biomass is standardized to kg per meter substrate length.

TABLE 1

Table 1. Total Nitrogen yield (tons) per hectare from the 2017 and 2018 growth cycle by potential harvest dates of belt spat collector treatments and alternative technologies.

FIGURE 4

Figure 4. (A) Mean and standard deviation (bars) of mussel biomass per meter length substrate for belt spat collector treatments over the growing season. Note that ice cover was present between late February and late March 2018 (gray shading). Farms are represented by color; SKIV (blue), DV (red), SALL (green). (B) Mean and standard deviation (bars) of mussel biomass per meter length substrate for alternative technologies and one belt spat collector treatment over the growing season in 2018. Farms are represented by color; SKIV (blue), DV (red), SALL (green).

Harvests

At harvest, total biomass was measured by treatment and line/tube in conjunction with sampling to evaluate the accuracy of the method. In 2017, sampling efforts sufficiently approximated biomass on spat collectors at both farms compared to harvest weights (DV, R² = 0.95, slope = 1.098, int = 0.98) (SKIV, R² = 0.80; slope = 0.932, int = 1.61). In 2018, spat collector biomass estimates by sampling approximated harvest weights, however, estimates on nets were complicated by harvest and cultivation issues, but corresponded well to measured harvested material.

Model Farm Production and Nutrient Extraction

Average winter harvest biomass yields and total N extraction per hectare, by technology and farm, for the two growth cycles are presented in Figure 5. Yields per hectare, by technology and farm are presented in Table 1 with statistics. Fractions of N and P in tissues, shells, and byssus are presented in Supplementary Table A3 (2017 growth cycle) and Supplementary Table A4 (2018 growth cycle), along with their respective proportions of total biomass at harvest. Potential extractive masses of each N and P are scaled to the model farm yield for each treatment. In general, total extractive potential of each substrate density and length configuration in December harvests were similar or greater than early-spring harvests (Supplementary Table A3). This trend is partially reflected in the biomass observations above and marginal changes or losses of N-content and/or tissue mass after the winter (in SKIV and SALL); nutrient content and dry weight yields are presented in Supplementary Table A5. While the conjunction of tissue yield, relative nutrient content in tissues, and total mussel biomass yields drive total nutrient extractive potential, biomass yield alone will account for >90% of the variation in total extractive potential. The critical metric in biomass yield is essentially total tissue yield, as tissue N, for example, accounts for 70.3 ± 2.5% of total N (2017 data, Supplementary Table A3). Losses of total biomass through the winter predominantly accounted for loss of nutrient extractive potential. Winter (November/December) harvests of 30 2 m belt spat collectors generally yielded ∼88 t mussels ha^–1 (1240 kg N ha^–1), 17.5 mesh net yielded ∼180 t mussels ha^–1 (2430 kg N ha^–1), and 25 cm mesh net yielded 154 t mussels ha^–1 (2070 kg N ha^–1) (Supplementary Table A3) between all sites. Including all biomass material (tissue, shells, and byssus), roughly 13.7 kg N and 0.9 kg P t^–1 harvested biomass can be extracted from the environment.

FIGURE 5

Figure 5. Upper panel: Biomass yield per model farm (t) by cultivation treatment. Lower panel: Nitrogen extraction (left axis) and phosphorus extraction (right axis, bars hashed) per hectare by cultivation treatment. Farms are represented by color; SKIV (blue), DV (red), SALL (green). Error bars represent standard deviations.

Belt Spat Collector Treatments

In terms of biomass yields, extractive potential, and buoyancy control, 30 2 m belt spat collectors provided optimal yields at the farm scale (Table 1 and Supplementary Table A3). In 2017, the 3 m belt spat collector treatments consistently exhibited significantly lower biomass yields and number of mussels per meter substrate than 2 m collectors (Supplementary Table A1). In shallow eutrophic estuaries, insufficient buoyancy can lead to spat collectors contacting the sea floor and exposing mussels to both predators and sulfide-rich sediments as was observed with the 3 m collectors. In September of 2017, the 3 m line at SALL broke and the majority of mussels were lost, hence no data for 3 m is presented in Table 1 or in the Appendix for SALL. Spacing treatments were similar in biomass yield throughout the season and during projected harvest periods. Mussel length and condition varied over the projected harvest periods, but typically followed a similar trend (Supplementary Table A2); 2 m collector treatments were similar and greater than both 3 m treatments. Biomass yields were highest in December relative to October and March; for 2 m collectors (Supplementary Table A3), these differences were significant only between the December and March periods (p < 0.05). December yields from the denser configuration (30 2 m) were 124% of March yields in SKIV, 103% in DV, and 113% in SALL (Supplementary Table A3). Belt spat collectors in 2018 demonstrated even higher yields than 2017 at similar size and condition (Supplementary Table A4). Accordingly, dense configuration of 2 m belt spat collectors, harvested in the winter, exhibit consistent high yields.

Alternative Technologies

In 2018, biomass yield followed substrate density, 17.5 cm nets providing the highest total yields overall (Table 1 and Supplementary Table A4). Mussel condition among the treatments followed an inverse trend to biomass, mussels on 17.5 cm nets were consistently smaller and had reduced meat content relative to other treatments (Supplementary Table A2, 2018). However, provided increased surface area and subsequent number of mussels on the surface, total biomass and nutrient yields were generally higher with increasing total substrate (Table 1 and Supplementary Table A4). In comparison to belt spat collectors, nets yielded considerably higher potential total biomass, 17.5 cm nets yielding 138% in SKIV, 199% in DV, 228% in SALL (Table 1 and Supplementary Table A4). Wider mesh nets (25 cm) also outperformed belt spat collectors, with 113% in SKIV, 142% in DV, and 229% in SALL. The two additional treatments in SALL, ladders and 20 cm net, yielded 131 and 258% of belt spat collectors, respectively; however, treatment replication was low. Comparably lower yields on nets in SKIV in 2018 to the other sites is mainly attributed to two storms which most likely caused sloughing of mussel aggregates, since large sections missing previously attached mussels was observed. Additionally, lower yields may be attributed to greater surface area and subsequently number of mussels within the specific section of the farms, as specific sections within each farm were utilized for both growth cycles (ruling out inter-farm locational preferentiality in growth potential). Nevertheless, the highest degree of replication was carried out at SKIV, and as such, conditions representative of a full-scale mitigation farm. Over all treatments, harvest periods (winter/spring), and farms, in a mixed model analysis (R² = 0.69, RMSE = 0.35), all treatments (p < 0.001) and harvest periods (p = 0.006) had significant effects on yield per hectare, while farm sites as random effects were non-significant (p > 0.05).

Substrate Normalized Yields

A consistent trend was observed for the two growth seasons and among different configurations of substrates following a proportional increase of biomass with substrate quantity, with diminishing gains approaching an asymptote. As materials between technologies have different surface area, total surface area per meter of longline or tube-net is regressed against biomass yield, Figure 6 (left panel) (R² = 0.48, p < 0.001), and mean shell length, Figure 6 (right panel) (R² = 0.32, p < 0.001). All cases exhibit a stronger logistic relationship with biomass yield than linear regression (R² = 0.32), indicating yields are positively impacted by substrate quantity up to a certain point, thereafter, increasing substrate quantity provides reduced gains in yield. In terms of total surface area per line/net, the first derivative of the non-linear function is less than 1 at approximately 6.1 m² substrate per meter line. This suggests only modest gains can be made by further increasing substrate coverage from those tested here. Shell length exhibits an inverse relationship with yield (Supplementary Table A1 and Figure 6), suggesting density-dependent growth in these tested environmental conditions. While correlations are relatively weak, these trends preliminarily suggest a relationship between growth, substrate area, and number of individuals in a given space; which will be regulated by self-thinning.

FIGURE 6

Figure 6. Left panel: Biomass yield per m longline or tube (kg) regressed against substrate surface area (m²) per m longline or tube at winter harvest [y = 309.33–451.53e^(-1.005x)], (R² = 0.48, p < 0.001). Right panel: Mean mussel shell length (mm) regressed against substrate surface area per meter line or tube (m²) at winter harvest [y = 45.46–1.09^∗e^(0.627)], (R² = 0.34, p < 0.001).

Spatial Analyses

Distributions of biomass and Condition Index at harvest in both DV and SKIV were plotted by sampling position and analyzed for directionality and clustering. Spatial autocorrelation (Moran’s Index) values and their p-values for each farm are presented in Table 2. DV was not analyzed in 2018 due to the low number of replicate lines and tube-nets. Distributions of Condition Index in SKIV (2017, 2018) and DV (2017) are presented in Figure 7. Reduced biomass yields in the center of units was detected to a limited degree, treatment differentiation provided higher clustering associations. Mussel Condition Index was clustered corresponding to edges (high values) and centrality (low values), and while significant, spatial gradients were not well defined (Figure 7 upper left panel, right panel) or limited to minor extents of the farm (Figure 7 lower left panel). Increasing mussel density (alternative technologies) tended to enhance spatial patterns in regard to yield and condition (Table 2).

TABLE 2

Table 2. Spatial autocorrelation (Moran’s Index) of biomass per meter substrate (kg m^–1) and Condition Index distributed over the farm.

FIGURE 7

Figure 7. Map of distribution of Condition Index = DW (mg)/SL³ (cm), for all treatments over the experimental section of each farm (SKIV and DV). The top left panel represents SKIV in the 2017 growth cycle (belt spat collectors), the bottom left represents SKIV in the 2018 growth cycle (tube-nets), the right panel represents DV in the 2017 growth cycle (belt spat collectors). Black lines represent the longlines or tube-nets as oriented in each farm. Interpolation over the farm was conducted by ordinary kriging (Haas and Wackernagel, 2006).

Specifically, in 2017, both farms exhibited negative unidirectional quadratic polynomial regression gradients (p < 0.001) of mussel condition toward the center of each farm, with significant clustering in SKIV at the center and edges (p < 0.01, Figure 7 upper left panel) and in DV only at the northern and southern extents (p < 0.05, Figure 7 right panel). Clustering and gradients of biomass yields was not extensive or exhibiting strong directionality. In 2018, mussel condition was negatively associated with the center and northwestern portions of the farm with significant unidirectional quadratic polynomial regression gradients (p < 0.05) and localized clustering of high (southeast and northeast, p < 0.04) and low values (center, p < 0.05) (Figure 7, lower left panel). In general, biomass yield clustering was limited.

Discussion

Production Optimization

Utilizing conventional long lines and a high-density belt spat collector configuration, yields exceeding 1600 t per model farm (∼85 t ha^–1) can be realized in a 6-month growth cycle. Alternative technologies exhibit potential to increase yields; exceeding 2500 t (∼133 t ha^–1) on average (Supplementary Table A4). The only previously published study and documentation on mitigation culture was conducted at the same location as SKIV in the present study in 2010–2011 and reported 1100 t after 12 months (Petersen et al., 2014; Nielsen et al., 2016). Furthermore, morphometric and linear yield were similar to 2017 results for lower density substrate spacing (60 cm) in the present study (October/November: 855 t vs. 1061 t in the present study). Higher potential yields were demonstrated here by increasing substrate density with conventional spat collectors (164% in November, 191% in December). When controlling for collector depth, high density of 2 m belt spat collectors in Limfjorden’s environment does not yet seem to approach production carrying capacity, even at sites with relatively lower average concentrations of chl-a (DV and SALL). Further increases in areal efficiency of nutrient extractive capacities was observed with the alternative technologies; up to 302% biomass yield with 17.5 cm nets in comparison with the prior study. While 20 cm nets provided the highest potential yields at SALL, replication of these tests at other sites is required before definitively concluding on mesh size optimization. These findings suggest substrate configuration will be very impactful on estimates of nutrient extractive potential for water management programs. Despite higher food quantities in SKIV, yields between sites were comparable. This indicates that the eutrophic conditions in the greater Limfjorden system are suitable for expanded mitigation cultivation. Nevertheless, the highest replication of treatments in the alternative technologies (tube-nets) was conducted at SKIV, so it is expected that the potential yields reported at SKIV will drive further parameterization.

Nutrient Extractive Capacities

Previously reported extraction potentials for a commercial scale mussel mitigation farm in SKIV was 0.6–0.9 t N ha^–1 and 0.03–0.04 t P ha^–1 (Petersen et al., 2014). Employing similar farming techniques and technologies, in 2017 the SKIV site exhibited potentials at December harvest of 0.8–1.4 t N ha^–1 and 0.06–0.1 t P ha^–1. As this production method in this environment aims at a single harvest per year, these extraction potentials are equivalent to an annual rate of extraction. Extractive potentials at DV and SALL were similar in magnitude, indicating common extractive capacities for Limfjorden waters employing marginally modified techniques from those that are currently utilized for human consumption. The alternative techniques provided higher extractive capacities, proportional to higher biomass yields. Tube-net systems, from all farms tested here, more than double potential extraction to 1.63–3.01 t N ha^–1 and 0.1–0.17 t P ha^–1. Nitrogen extraction capacities at this magnitude are similar to efficient constructed wetlands (Vymazal, 2007) and integrated buffer zones (Zak et al., 2018), while P removal was generally lower.

Farm Substrate Configuration and Management

From experiences with spacing and length experiments in 2017, this study suggests that denser spacing is more advantageous than lengthening the collectors to take advantage of an increased proportion of the water column to increase areal efficiency. In shallow estuaries, there is reduced buffer distance between the bottom of the suspended canopy and the benthic sediment surface. The increased weight of mussels on longer collectors requires proportionally greater buoyancy maintenance as loads on the mainline increase more rapidly as mussels grow. Buoyancy is added in moderation, as excessive buoyancy in conjunction with high frequency wind-driven wavelets or small waves can lead to displacement and loss of mussels (Young, 1985; Carrington et al., 2015). As farmers were typically maintaining buoyancy according to a schedule corresponding to the conventional configuration (i.e., 2 m depth and 60 cm spacing), the 3 m collectors in this study consistently sank below stable buoyancy depth and were occasionally in contact with the sea floor (water depths of 5.5–7 m). Contact with the sea floor will exposes the mussels to predators e.g., Asterias rubens, Carcinus maenas (Kamermans et al., 2009) and sulfide-rich organic sediments, resulting in reduced overall mussel yield and lower mussel condition. This is evident in lower than expected yields from 3 m collectors in this study, and such predation and degradation conditions were observed on these collectors. As such, buoyancy strategies could be employed to sufficiently maintain longer collectors, which should increase potential yields from those documented in this study but would increase operational costs relative to shorter collectors.

Further increasing substrate area requires modifications in buoyancy application. The ladder system requires reduced buoyancy maintenance as total buoyancy is allocated across fewer individual buoys. Higher biomass yields per line were observed with ladders relative to belt spat collectors at SALL, suggesting this technology may exhibit higher nutrient extractive potential given the same area, and follows the surface area: yield trend, however, should be replicated in further study. Establishment of mitigation farms in deeper waters could find advantage in the adoption of these systems.

Tube-net configurations (i.e., Smartfarm system) require little to no buoyancy maintenance over the growth period, unless mussels are expected to be overwintered, where tube-nets provide the least flexibility in buoyancy control. Protecting these units from ice coverage in shallow areas (such as experienced in February-March 2018) requires exchanging air for water within the tube, resulting in the unit – and attached biota – sinking to the sea floor, which will lead to loss of biomass to soft-sediment exposure and/or predation (e.g., Kamermans et al., 2009). Harvest before potential ice coverage forgoes these issues. In comparison with both conventional belt spat collectors at high density (30 cm) and ladders, tube-nets provided considerably higher potential yields at early winter harvest; at best exceeding 300%. The capital investments required for tube-nets are much higher than long line systems (∼190%), however, only require ∼10% of the maintenance labor costs (Filippelli et al., in review). The much higher potential yields may justify higher initial costs, however, this requires further investigation as equipment and mitigation mussel market forces will likely adjust according to the scale of industrial development, potential valorization of ecosystem services, and future trends in feeds markets (Asche et al., 2013).

Timing

Harvest dates and sampling indicated, in general, that December harvests provided the highest potential biomass yields and nutrient extractive potential relative to early spring of the following year. This trend is reflected in the prior study assessing production capabilities of a mitigation farm, where October yields were 136% of March yields (Nielsen et al., 2016). Ice cover in that study period persisted for three months, while in the present study, ice cover lasted only 1 month. May harvests were shown to exhibit the highest potential yields in the previous study, while operational requirements for farm management mandated harvest before May in the present study. In the present study, overwintering mitigation mussels for an early spring harvest provided little advantage in contrast to a November or December harvest. With added costs to maintain the units over the winter months, loss of biomass, and subsequent lack of somatic accretion due to harvest prior to spring phytoplankton blooms, harvests initiated in the later autumn and continuing through the winter are favorable in terms of mitigating nutrient loads. Operationally, an extended harvest over the winter is suitable provided ice cover or increasingly harsh winter weather does not compromise retention of biomass, and provide additional time to maintain farm materials and equipment prior to the next settlement in late-May to June. As the capacity for processing large volumes of materials is currently limited, an extended harvest season would be required with the proliferation of farms in the Limfjorden. Expansion of mitigation production is therefore likely to require an accompanied expansion of processing capabilities and infrastructure.

Furthermore, it is also important to consider seasonal and related variability in nutrient dynamics (i.e., cycles and loads) when maximizing extractive potential. Accumulated sedimentation, mineralization, and higher ammonium excretion in the spring may reduce the net extractive capacity of mitigation units when taking an ecosystem budgetary approach (Holmer et al., 2015). As such, and with little extractive gain by overwintering, from an integrated production and ecosystem perspective, 6–8 month harvest cycles are optimal for mitigation cultivation in the present environmental contexts.

Variability, Local Effects, and Limitations to Production

Prior work has demonstrated an inverse relationship between mussel growth and density within a mussel farm (Cubillo et al., 2012). As increased surface area increases total population in a given space, integrated filtration pressure can lead to food limitation at the boundary layer around the mussel aggregates or at a larger scale further afield (Petersen et al., 2019b). Such a pattern emerges when comparing condition or shell length between treatments (Supplementary Table A2). In both SKIV and DV, shell length and condition differ significantly between spat collectors and the finer mesh net (17.5 cm) (Supplementary Table A2, 2018 section).

It can be assumed that maximizing the substrate in a given mitigation farm will simply provide the greatest nutrient extractive potential. This is represented in the case of the longline and belt spat collector treatments, but when normalizing total yield per longline or net to surface area (Figure 6), a diminishing return on yields is observed to an asymptote, namely, in between the 17.5 and 25 cm nets. This is due largely to two factors: (1) at higher substrate quantities, there is less space for mussel aggregate expansion, and (2) overlapping food depletion fields reduce interspecific growth (Saurel et al., 2013; Petersen et al., 2019b). Reduced growth will be increasingly pronounced in estuaries with lower mean food concentrations, lower flux of food, or if cumulative depletion of basin-scale food resources induces similar conditions (Strohmeier et al., 2008).

Spatial patterns can expose trends related to the interaction of the farm structure with the local biophysical regime and interaction among treatments within the farm. Spatial clustering of high and low biomass and mussel condition indicates farm-scale localized food depletion and is likely to be detected at downstream or central sections, depending on currents and food concentrations. Spatial distribution of biomass yields was weakly oriented in both farms, indicating food flux was sufficient for overall production. While reduced condition at harvest in the center of the units in both SKIV and DV was observed (Table 2 and Figure 7), these spatial trends were detected between treatments, as well as independent of treatment. This implies both sufficient spatial randomization of treatment assignment, and the spatial orientation of reduced condition is a feature of integrated depletion effects. Mussel condition is directly related to overall extractive potential as the majority of nutrients are bound in the soft tissues provided differences in size (shell length) and numbers (density) are negligible; directional dependency of condition can inform on modified harvest strategies or adjustments to the physical configuration. Greater mussel condition was observed in larger mesh nets (Supplementary Table A2) and can be likely attributed to reduced micro-scale depletion fields as mesh holes were not entirely obstructed by mussels, permitting food flux through the net structure. It has been demonstrated that chl-a concentrations decline rapidly within a region of 5–10 cm from the mussel aggregate surface (Petersen et al., 2019b), which integrates into a larger volume for smaller mesh nets. When considering predominant current directions and integrated depletion effects at both the aggregate surface and farm-scale, clustering magnitudes (low and high values) and ranges were positively associated with substrate density (Table 2).

The critical metric in mitigation culture is total nutrient extractive potential. This manifests as total summary biomass, which is largely commensurate to total tissue content. In both of the commercial farms (SKIV and DV), total biomass yield, total N, and total P were greatest in the case of the 17.5 cm nets (Figure 5). This indicates that at the farm-scale, depletion-induced concavity in total yield or nutrient extraction had not exceeded production carrying capacity. Nonetheless, expansion of production will require higher resolution analysis linking micro-scale depletion dynamics with basin-scale biophysical dynamics influencing depletion. Additionally, it cannot be understated that farm management, specifically in relation to maintenance of biomass retention, will markedly influence yields.

Future Considerations

If mitigation farms proliferate, compounding large-scale effects on phytoplankton is expected as an outcome in eutrophic waters. Basin-scale reduction in food supplies and quality will influence farm configuration, as the flux of food in a reduced gradient will support less biomass in upstream units if situated in close proximity. Furthermore, it has been demonstrated in heavily cultivated waters that ecological space is broadened for smaller phytoplankton classes (Froján et al., 2018; Cranford, 2019); however, in eutrophic or some extensively cultivated basins, smaller plankton classes may be increasingly utilized under field conditions (Sonier et al., 2016) or transported through the food chain to larger classes (Haraguchi et al., 2018). As of 2019, 7.34 km² of the Limfjorden are permitted for cultivation of mussels, where the overall surface area 1526 km² (Dinesen et al., 2011), or 0.5% of the fjord system; which suggests room for expansion in this highly productive system, in comparison with more extensively utilized areas (e.g., St. Peter’s Bay, Canada, 38% utilized; Sonier et al., 2016). Nevertheless, ecological impacts of expanded and intensive mitigation cultivation require further investigation through field study and modeling efforts to estimate impacts of increased scale. Basin-scale configuration of farms to optimize nutrient extraction will require strategic placement that balances (1) biophysical forcing, (2) multiple-use demands, (3) existing abatement mechanisms in the watershed, (4) ecological sensitivity of the benthos, (5) local capacity in farm technical management, and (6) consideration of impacts on natural recruitment in mussel beds (Molinet et al., 2017).

As mussels immobilize organic matter from the water column as fecal and pseudofecal matter, there is potential for organic enrichment of sediments directly within and around the farm. Sedimentation rates of organic matter have been demonstrated to be higher within a farm than at reference positions within the same water body (Carlsson et al., 2009), however, basin-scale sedimentation is reduced in eutrophic water bodies (Timmermann et al., 2019). It is expected that further commercial scale, in situ investigation of environmental conditions and farm configuration will clarify questions of nutrient budgeting.

The fate of mitigation mussels post-harvest is currently undergoing technical³ and economic investigation (Filippelli et al., in review). While many mussels will be appropriately sized for the human consumption market, utilization in feeds or otherwise will follow market and technological conditions. Further elaboration on ecosystem services, ecological impacts, and economics of production of mitigation mussels are discussed in reviews by Petersen et al. (2016, 2019a).

Conclusion

Increased areal efficiency of mitigation mussel production was demonstrated in this study by modifying conventional cultivation methods and adoption of alternative technologies. Nutrient extractive potential was increased over prior work and exhibits rates similar to highly efficient constructed terrestrial mitigation mechanisms. While sites in this study provided different conditions for growth, total yields were similar and suggest Limfjorden waters are appropriate for further expansion of mitigation culture. Further commercial-scale demonstration of alternative technologies would facilitate heightened precision and accuracy of harvest potentials under different environmental conditions. As mitigation cultivation is intended to augment other nutrient abatement mechanisms, and emphasis on targeted approaches is becoming standard practice, extraction capacity of a given farm or collection of farms within a sub-basin needs to be linked to local ecological conditions. The magnitude of filtration pressure and potential summary impacts on nutrient fates implies that mitigation farming will modify local and meso-scale conditions. Similar to other abatement tools, response to variability in conditions over the year, site-specific attributes, and cultivation practices will affect extraction efficiency. Further study is required to document farm and basin-scale impacts on phytoplankton and organic matter distribution relative to strategic deployment of mitigation farms.

Data Availability Statement

All datasets generated for this study are included in the article/Supplementary Material.

Author Contributions

DT, CS, PN, and JP contributed to conception and design of the study. DT and CS contributed to data curation. DT contributed to statistical and formal analyses, investigation, and writing of the original draft. CS, PN, and JP contributed to funding acquisition, project administration, resources, supervision, and review and editing of the manuscript. All authors contributed to the manuscript revision, read and approved the submitted version.

Funding

This study was conducted under the MuMiPro project (Mussel farming, mitigation, and protein source for organic husbandry) project, 6150-00008B, funded by the Danish Innovation Fund, and the BONUS OPTIMUS project (Optimization of mussel mitigation cultures for fish feed in the Baltic Sea) as part of the, funded by the BONUS program (Baltic Organisations’ Network for Funding Science EEIG, Art185), funded jointly by EU and Innovation Fund Denmark.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Finn Bak and Niels-Peter Nielsen for their dedication to the project work, technical expertise, and contributions to project execution. We further acknowledge the diligent assistance of Danish Shellfish Center staff in sampling and sample processing: Anita Hansen, Helge Boesen, Jacob Emil Jensen, Kasper Lenda Andersen, Lars Bach Jensen, Lars Kyed Andersen, and Lotte Shrøder. We thank Ivar Lund and Ulla Sproegel for nutrient analyses. We also thank the section editor and the two reviewers for their constructive feedback on the manuscript, improving the quality of the current version. Lastly, we acknowledge Seafood Limfjord and Wittrup Seafood for use of their farms, personnel, equipment, and cooperation in research.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2019.00698/full#supplementary-material

Footnotes

References

Andersen, H. E., Blicher-Mathiesen, G., Bechmann, M., Povilaitis, A., Iital, A., Lagzdins, A., et al. (2014). Mitigating diffuse nitrogen losses in the Nordic-Baltic countries. Agric. Ecosyst. Environ. 195, 53–60. doi: 10.1016/J.AGEE.2014.05.009