Modelling mussel (Mytilus spp.) microplastic accumulation

Microplastics (MPs) are a contaminant of growing concern due to their widespread distribution and interactions with marine species, such as filter feeders. To investigate the MPs accumulation in wild and cultured mussels, a dynamic energy budget (DEB) model was developed and validated with the available field data of Mytilus edulis (M. edulis, wild) from the North Sea and Mytilus galloprovincialis (M. galloprovincialis, cultured) from the northern Ionian Sea. Towards a generic DEB model, the site-specific model parameter, half-saturation coefficient (X_k), was applied as a power function of food density for the cultured mussel, while for the wild mussel it was calibrated to a constant value. The DEB-accumulation model simulated the uptake and excretion rate of MPs, taking into account environmental characteristics (temperature and chlorophyll a). An accumulation of MPs equal to 0.53 particles per individual (fresh tissue mass 1.9 g) and 0.91 particles per individual (fresh tissue mass 3.3 g) was simulated for the wild and cultured mussel after 4 and 1 years respectively, in agreement with the field data. The inverse experiments investigating the depuration time of the wild and cultured mussel in a clean-from-MPs environment showed a 90 % removal of MPs load after 2.5 and 12 d respectively. Furthermore, sensitivity tests on model parameters and forcing functions highlighted that besides MPs concentration, the accumulation is highly dependent on temperature and chlorophyll a of the surrounding environment. For this reason, an empirical equation was found, directly relating the environmental concentration of MPs, with the seawater temperature, chlorophyll a, and the mussel's soft tissue MPs load.

1 Introduction

Microplastic (MP) particles are synthetic organic polymers with size below 5 mm (Arthur et al., 2009) that originate from a variety of sources including the following: those particles that are manufactured for particular household or industrial activities, such as facial scrubs, toothpastes, and resin pellets used in the plastic industry (primary MPs), and those formed from the fragmentation of larger plastic items (secondary MPs) (GESAMP, 2015). Eriksen et al. (2014) estimated that more than 5 trillion microplastic particles, weighing over 250 000 t, float in the oceans. Due to their composition, density, and shape, MPs are highly persistent in the environment and are, therefore, accumulating at different rates in different marine compartments: surface and deeper layers in the water column, as well as at the seafloor and within the sediments (Moore et al., 2001; Lattin et al., 2004; Thompson, 2004; Lusher, 2015). Since the majority of MPs entering the marine environment originate from the land (i.e. landfills, littering of beaches and coastal areas, rivers, floodwaters, untreated municipal sewerage, and industrial emissions), the threat of MPs pollution in the coastal zone puts considerable pressure on the coastal ecosystems (Cole et al., 2011; Andrady, 2011). In recent years, initiatives under various projects (i.e. CLAIM, DeFishGear) aim at evaluating the threat and impact of marine litter pollution; the European framework of JERICO-RI focuses on a sustainable research infrastructure in the coastal area to support the monitoring, science, and management of coastal marine areas (http://www.jerico-ri.eu/, last access: 27 July 2020). In the framework of JERICO-NEXT, a recent study addressed the environmental threats and gaps with monitoring programmes in European coastal waters, including the marine litter (i.e. MPs), as one of the most commonly identified threats to the marine environment and highlighted the need for improved monitoring of the MPs distribution and their impacts in European coastal environments (Painting et al., 2020).

Numerous studies have revealed that MPs are ingested either directly or through lower trophic prey by animals at all levels of the food web – from zooplankton (Cole et al., 2013), small pelagic fish, and mussels (Digka et al., 2018a) to mesopelagic fish (Wieczorek et al., 2018) and large predators like tuna and swordfish (Romeo et al., 2015). Microplastic ingestion by marine animals can potentially affect animal health and raises toxicity concerns, since plastics can facilitate the transfer of chemical additives and/or hydrophobic organic contaminants to biota (Mato et al., 2001; Rios et al., 2007; Teuten et al., 2007, 2009; Hirai et al., 2011). Humans, as top predators, are also contaminated by MPs (Schwabl et al., 2019). Mussels and small fish that are commonly consumed whole, without removing digestive tracts, where MPs are concentrated, are among the most likely pathways for MPs to embed in the human diet (Smith et al., 2018). Especially regarding marine organisms (i.e. mussels), it is notable that the levels of their contamination have been added to the European database (http://www.ecsafeseafooddbase.eu, last access: 27 July 2020) as an environmental variable of growing concern, reflecting the health status (Marine Strategy Framework Directive, MSFD, Descriptor 10 – Marine Litter; Decision 2017/848/EU) (De Witte et al., 2014; Vandermeersch et al., 2015; Digka et al., 2018a). Today, a series of studies have noted the presence of MPs in mussels' tissue intended for human consumption (Van Cauwenberghe and Janssen, 2014; Mathalon and Hill, 2014; Li et al., 2016, 2018; Hantoro et al., 2019). For instance, in a recent study, Li et al. (2018) sampled mussels from coastal waters and supermarkets in the UK and estimated that a plate of 100 g mussels contains 70 MPs that will be ingested by the consumer. The presence of MPs in mussels has been also demonstrated during laboratory trials in their faeces, intestinal tract (Von Moos et al., 2012; Van Cauwenberghe et al., 2015; Wegner et al., 2012; Khan and Prezant, 2018), and circulatory system (Browne et al., 2008). Other laboratory studies showed several effects of microplastic ingestion in laboratory-exposed mussels, including histological changes, inflammatory responses, immunological alterations, lysosomal membrane destabilization, reduced filtering activity, neurotoxic effects, oxidative stress effects, increase in haemocyte mortality, dysplasia, genotoxicity, and transcriptional responses (reviewed by Li et al., 2019). However, the tested concentrations of MPs in laboratory experiments are frequently unrealistic, being several orders of magnitude higher (2 to 7 orders of magnitude) than the observed seawater concentrations (Van Cauwenberge et al., 2015; Lenz et al., 2016).

Mussels, through their extensive filtering activity, feed on planktonic organisms that have a similar size to MPs (Browne et al., 2007), and considering also their inability to select particles with high energy value (i.e. phytoplankton) during filtration (Vahl, 1972; Saraiva et al., 2011a), they are directly exposed to MPs contamination. Recent studies suggest a positive linear correlation between MPs concentration in mussels and surrounding waters (Capolupo et al., 2018; Qu et al., 2018; Li et al., 2019). The filtering activity of mussels, which directly affects the resulting MPs accumulation, is a complicated process that is controlled by other factors (food availability, temperature, tides etc.).

The purpose of the present work is to study the accumulation of MPs in mussels and reveal relations between the accumulated concentrations in mussels' soft parts and environmental features. In this context, an accumulation model was developed based on dynamic energy budget theory (DEB; Kooijman, 2000) and applied in two different regions, in two different modes of life (wild and cultivated): in the North Sea (Mytilus edulis (M. edulis), wild) and in the northern Ionian Sea (Mytilus galloprovincialis (M. galloprovincialis), cultivated). DEB theory provides all the necessary detail to model the feeding processes and aspects of the mussel metabolism, taking into account the impact of the environmental variability on the simulated individual. Apart from modelling the growth of bivalves (Rosland et al., 2009; Sarà et al., 2012; Thomas et al., 2011; Saraiva et al., 2012; Hatzonikolakis et al., 2017; Monaco and McQuaid, 2018), DEB models have been used to study other processes as well, such as bioaccumulation of PCBs (polychlorinated biphenyls) and POPs (persistent organic compounds) (Zaldivar, 2008), trace metals (Casas and Bacher, 2006), and the impact of climate change on individual's physiology (Sarà et al., 2014). However, to our knowledge this is the first time that a DEB-based model is used to assess the uptake and excretion rates of MPs in mussels.

2 Materials and methods

2.1 Study areas and field data

The North Sea is a marginal sea on the continental shelf of north-western Europe with a total surface area of 850 000 km² and is bounded by the coastlines of nine countries. The sea is shallow (mean depth 90 m), getting deeper towards the north (up to 725 m), and the semidiurnal tide (tidal range 0–5 m) is the dominant feature of the region (Otto et al., 1990). Major rivers, such as the Rhine, Elbe, Weser, Ems, and Thames discharge into the southern part of the sea (Lacroix et al., 2004), making this area a productive ecosystem. In this study, the area is limited to along the French, Belgian, and Dutch North Sea coast (50.98–51.46^∘ N, 1.75–3.54^∘ W). This is located close to harbours, where shipping, industrial, and agricultural activity is high, putting considerable pressure on the ecological systems of the region (Van Cauwenberghe et al., 2015).

The MPs concentration in mussels' tissue and seawater that were used to validate and force the model respectively at its North Sea implementation were derived from Van Cauwenberghe et al. (2015). Van Cauwenberghe et al. (2015) examined the presence of MPs in wild mussels (M. edulis) and thus collected both biota and water at six sampling stations along the French, Belgian, and Dutch North Sea coast in late summer of 2011. M. edulis (mean shell length: 4±0.5 cm and wet weight (w.w.): 2±0.7 g) and water samples were randomly collected in the local breakwaters, in order to assess the MPs concentration in the organisms and their habitat. MPs were present in all analysed samples, both organisms and water. Seawater samples (N=12) had MPs (<1 mm) on average 0.4±0.3 particles L⁻¹ (range: 0.0–0.8 particles L⁻¹), and M. edulis contained on average 0.2±0.3 particles g⁻¹ w.w. (or 0.4±0.3 particles per individual) (Van Cauwenberghe et al., 2015). The size range of MPs found within the mussels was 20–90 µm.

The northern Ionian Sea is located in the transition zone between the Adriatic and Ionian Sea. The long and complex coastline presents a high diversity of hydrodynamic and sedimentary features. Rivers discharging into the northern Ionian Sea include Kalamas/Thyamis (Greece) and Butrinto (Albania) (Skoulikidis et al., 2009; Vlachogianni et al., 2017), making the area suitable for aquaculture. Small farming sites and shellfish grounds are operating in Thesprotia (northwestern Ionian Sea) (Theodorou et al., 2011). The main source of marine litter inputs in the area originates from anthropogenic activities that mainly include shoreline tourism and recreational activities, poor wastewater management, agriculture, fisheries, aquacultures, and shipping (Vlachogianni et al., 2017; Digka et al., 2018a). According to Politikos et al. (2020), the area around the Corfu island (northern Ionian Sea) is characterized as a retention area of litter particles probably due to the prevailing weak coastal circulation. Furthermore, a northward current on the east Ionian Sea facilitates the transfer of litter particles towards the Adriatic Sea, which has been characterized as a hotspot of marine litter and one of the most affected areas in the Mediterranean Sea (Pasquini et al., 2016; Vlachogianni et al., 2017; Liubartseva et al., 2018; Politikos et al., 2020).

The field data used to validate the model output in the N Ionian Sea were obtained from Digka et al. (2018a, b). In the framework of the “DeFishGear” project, mussels (M. galloprovincialis) were collected by hand from a longline-type mussel culture farm in Thesprotia (39.606567^∘ N, 20.149421^∘ E), in summer 2015 (end of June) at a sampling depth up to 3 m (Digka et al., 2018a). The average MPs accumulation was calculated from a total population of 40 mussels originated from the farm, with 18 of them were found to be contaminated with MPs (46.25 %). The average load of MPs (size <1 mm) per mussel (mean shell length 5.0±0.3 cm) was 0.9±0.2 particles per individual, and the size of MPs found in the mussel's tissue ranged from 55 to 620 µm. Both clean and contaminated mussels were included in the calculated mean value in order to represent the mean state of the contamination level for the individual inhabiting the study area.

The seawater concentration of MPs for the N Ionian Sea implementation was obtained from Digka et al. (2018b) and the DeFishGear project results (http://www.defishgear.net/project/main-lines-of-activities, last access: 27 July 2020). In total, 12 manta net tows were conducted in the region, collecting a total number of n₁=2027 particles on October 2014 and n₂=1332 on April 2015, leading to an average of 280 particles per tow with size <1 mm and >330 µm (Digka et al., 2018b). In order to estimate the mean MPs concentration in the region, expressed as particles per volume, the dimensions of the manta net (W: 60 cm, H: 24 cm, rectangular frame opening; mesh size 330 µm) and the sampling distance of each tow (∼2 km) were used by multiplying the sample surface of the net by the trawled distance in metres (Maes et al., 2017), which resulted in a mean MPs concentration of 1.17 particles m⁻³ (233 333 particles km⁻²). Moreover, in the wider region of the Adriatic Sea, Zeri et al. (2018) found a mean density of 315 009±568 578 particles km⁻² (1.58±2.84 particles m⁻³), out of which 34 % were sized <1 mm. A relatively high value of standard deviation (1 order of magnitude higher than the mean value) is adopted (0.0012±0.024 particles L⁻¹), considering that the mussel farm is established in an enclosed gulf and close to the coast, since, according to Zeri et al. (2018), the abundance of MPs is 1 order of magnitude higher in inshore (<4 km) compared to offshore waters (>4 km). Furthermore, it may be assumed that the adopted range (standard deviation is also multiplied by a factor of 2) includes also the smaller particles sized between 50 and <330 µm, which have been found in mussel tissue (Digka et al., 2018a) but were overlooked during the seawater sampling due to the manta net's mesh size (>330 µm). According to Enders et al. (2015) the relative abundance of small particles (50–300 µm) compared to particles larger than 300 µm is approximately 50 %.

2.2 DEB model description

In the present study, a DEB (Kooijman, 2000, 2010) model is used as basis to simulate the accumulation of MPs by mussels. In DEB theory (Kooijman, 2000), the energy assimilated through food by the simulated individual is stored in a reserve compartment from where a fixed energy fraction κ is allocated for growth and somatic maintenance, with a priority for maintenance. The remaining energy (1−κ) is spent on maturity maintenance and reproduction. The individual's condition is defined by the dynamics of three state variables: energy reserves E (joules), structural volume V (cm³), and energy allocated to reproduction R (joules). The energy flow through the organism is controlled by the fluctuations of the available food density and temperature characterizing the surrounding environment.

The DEB model implemented here is an extended version of the model described in Hatzonikolakis et al. (2017), where the growth of the Mediterranean mussel is simulated taking into account only the assimilation rate of the individual. Since the present study focuses on the MPs accumulation, it is crucial to include a detailed representation of the mussel's feeding mechanism. In this context, the DEB model was extended by including the clearance (C_r), filtration (${\stackrel{\mathrm{´}}{p}}_{\mathrm{XiF}}$), and ingestion (${\stackrel{\mathrm{´}}{p}}_{\mathrm{XiI}}$) rates of the mussel, following Saraiva et al. (2011a), assuming that all parameters referred to as silt (or inedible particles) are applicable also to MP particles. In this approach, a pre-ingestive selection occurs between filtration and ingestion, returning the rejected material in the water through pseudofaeces (J_pfi). Consequently, energy is assimilated through food, and the non-assimilated particles are excreted through the faeces production (J_f). The model's equations, variables, and parameters are shown in Tables 1, 2, and 3 respectively. The scaled functional response f (Eq. 5, Table 1), which regulates the assimilation rate, is modified following Kooijman (2006) to include an inorganic term representing the non-digestible matter, i.e. microplastics: $f=X/\left(X+K_y\right)$ and $K_y=X_K\cdot (\mathrm{1}+Y/Y_K$), where Y and Y_k are the concentration of MPs, converted from particles per litre to grams per cubic metre (Everaert et al., 2018) and the half-saturation coefficient of inorganic particles here represented by MPs (g m⁻³) respectively. Thus, the assimilation rate that is regulated by f is decreasing when the concentration of MPs is increased. The same approach is followed by other authors who considered inedible particles in the mussel's diet (Ren, 2009; Troost et al., 2010). During the filtration process the same clearance rate for all particles is used $(\left\{{\dot{C}}_{\mathrm{R}}\right\}$), representing the same searching rate for food that depends on the organism maximum capacity $(\left\{{\dot{C}}_{\mathrm{Rm}}\right\}$) and environmental particle concentrations (Vahl, 1972; Widdows et al., 1979; Cucci et al., 1989). During the ingestion process the mussel is able to selectively ingest food particles and reject inedible material, in order to increase the organic content of the ingested material (Kiørboe and Møhlenberg, 1981; Jørgensen et al., 1990; Prins et al., 1991; Maire et al., 2007; Ren, 2009; Saraiva et al., 2011a). This selection is reflected by the different binding probabilities adopted for each type of particle (ρ₁ for algae particles and ρ₂ for inorganic particles, i.e. MPs; see Eq. 14 and Table 3). The equations representing the feeding processes handle each type of particle separately, while there is interference between the simultaneous handling of different particle types (Eqs. 12–14, Table 1) (Saraiva et al., 2011a). Finally, during the assimilation process, suspended matter (i.e. MPs) that the mussel is not able to assimilate due to its different chemical composition from the reserve compartment (Saraiva et al., 2011a) or incipient saturation at high algal concentrations (Riisgård et al., 2011) results in the faeces production (Eq. 16, Table 1).

Table 1 Dynamic energy budget model: equations. See Table 2 for model variables, Table 3 for parameters, and Table 5 for initial values.

^* Notation refers to feeding equations handling each type of suspended matter separately (i=1 for algae and i=2 for microplastics) where unit transformation is applied when it is necessary (see Table 3).

2.3 Microplastics accumulation submodel

With the DEB model as a basis, a submodel describing the MPs accumulation by the mussel was developed, assuming that the presence of MPs in the ambient water does not cause a significant adverse effect on the organisms' overall energy budget, in accordance with laboratory experiments, conducted in mussel species (Van Cauwenberghe et al., 2015, for Mytilus edulis; Santana et al., 2018, for the mussel Perna perna). Additionally, it was assumed that the mussel filtrates MPs present in the water, without the ability of selecting between the high-energy-valued particles and the MPs during the filtration process (Van Cauwenberghe et al., 2015; Von Moos et al., 2012; Browne et al., 2008; Digka et al., 2018a among others). The uptake of MPs from the environment is taken into account through the process of clearance/filtration rate, while the excretion of the contaminant is derived from two processes: (i) pseudofaeces production and (ii) faeces production. The resulting MPs accumulation is influenced by external environmental factors (MPs concentration, food availability, temperature) and internal biological processes (clearance, filtration, ingestion, growth). The following differential equation describes the change of the individual MPs accumulation (C, particles per individual), taking into account the processes mentioned above:

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{\mathrm{d}C}{\mathrm{d}t}}=C_{\mathrm{env}}\cdot {\stackrel{\mathrm{´}}{C}}_{\mathrm{R}}-{\stackrel{\mathrm{´}}{J}}_{\mathrm{pf}\mathrm{2}}-k_{\mathrm{f}}\cdot {\displaystyle \frac{{\stackrel{\mathrm{´}}{J}}_{\mathrm{f}}}{p_{X_{\mathrm{1}}I}}}\cdot C,\end{array}$$

where ${\stackrel{\mathrm{´}}{C}}_{\mathrm{R}}$ is the clearance rate for water (L h⁻¹), containing a concentration of MPs C_env (particles L⁻¹). The terms of ${\stackrel{\mathrm{´}}{J}}_{\mathrm{pf}\mathrm{2}}$ and $\frac{{\stackrel{\mathrm{´}}{J}}_{\mathrm{f}}}{p_{X_{\mathrm{1}}I}}$ represent the elimination rate of MPs through pseudofaeces (particles h⁻¹) and the nondimensional rate of faeces production with respect to the ingestion rate respectively (see Table 1, Eqs. 15–16). The parameter k_f represents the post-ingestive selection mechanism utilized by the mussel to incorporate indigestible material (i.e. MPs) into faeces and was calibrated using the available field data of mussel MPs accumulation from both study areas (Table 3). A mussel is able to discriminate among particles in the gut based on size, density and chemical properties of the particles (i.e. between microalgae and inorganic material) and thus to eliminate them through faeces (Ward et al., 2019a, and references therein). In this context, the pseudofaeces production incorporates the rejected MPs prior to the ingestion, while the faeces production includes MPs that are rejected along with the food particles that are not assimilated by the mussel. The model's time step has been set to 1 h in order to capture the dynamics of the rapidly changing processes, such as feeding and excretion.

Table 2 Dynamic energy budget model: variables

Table 3 Dynamic energy budget model: parameters.

^* Units of moles of carbon (mol C) converted to milligrams of CHL a (mg CHL a) by multiplying with the factor $\frac{\mathrm{12}\times {\mathrm{10}}^{\mathrm{3}}}{\mathrm{50}}$ assuming carbon:CHL a ratio of 50 (Hatzonikolakis et al., 2017).

2.4 Environmental drivers

Besides MPs concentration in the seawater, the DEB model is forced by sea surface temperature (SST) and food availability, represented by chlorophyll a concentrations (CHL a, an index of phytoplankton biomass). M. edulis has been demonstrated to filter suspended particles greater than 1 µm, a size class that includes all of the phytoplankton, zooplankton, and much of the detritus (Vahl, 1972; Møhlenberg and Riisgård, 1978; Saraiva et al., 2011a; Strohmeier et al., 2012), including even aggregated picoplankton-size particles (i.e. marine snow) (Kach and Ward, 2007; Ward and Kach, 2009). CHL a has been considered the most reliable food quantifier for the calculation of DEB shellfish parameters (Pouvreau et al., 2006; Sarà et al., 2012; Hatzonikolakis et al., 2017, and references therein). Hatzonikolakis et al. (2017) have tested the performance of the model, considering also particulate organic carbon (POC) in the mussel's diet, which, however, did not have an important impact on the model's skill against field data in the Mediterranean Sea study areas. This outcome agrees with Troost et al. (2010) demonstration that POC contributes to the mussel's diet when CHL a concentrations are low in the southwest of Netherlands. Thus, in the present study, only CHL a is considered as the available food source for mussels originated from the southern North Sea and the northern Ionian Sea. For both study areas SST and CHL a are derived from daily satellite data, a method also used by other authors (i.e. Thomas et al., 2011; Monaco and McQuaid, 2018).

In the North Sea, SST data were obtained from daily satellite images provided by Copernicus Marine Environmental Monitoring Service (CMEMS) at 0.04^∘ spatial resolution. CHL a data obtained from the Globcolour daily multi-sensor product provided by CMEMS at 1 km spatial resolution, based on the OC5 algorithm of Gohin et al. (2002) (http://marine.copernicus.eu/, last access: 27 July 2020, generated using CMEMS Products, production centre ACRI-ST). The environmental forcing data (SST, CHL a) were averaged over the study area (51.08–51.44^∘ N, 2.19–3.45^∘ E), covering the period 2007–2011 (5 years), in order to realistically simulate the wild mussel's growth harvested in late summer 2011 (Van Cauwenberghe et al., 2015). It is notable that the study area of the North Sea belongs to Case II waters (coastal region), where algorithms tend to overestimate CHL a concentrations. In optically complex Case II waters, CHL a cannot readily be distinguished from particulate matter and/or yellow substances (dissolved organic matter), and so global chlorophyll algorithms are less reliable (IOCCG, 2000). However, the CHL a dataset that was used was found in good agreement with available in situ data from the ICES database (https://www.ices.dk/data/Pages/default.aspx, last access: 27 July 2020) for the specific study area and time period (Fig. 1), showing a relatively smaller bias and better time–space coverage, as compared with other tested remote sensing datasets (not shown) (i.e. regional chlorophyll product available for the North West Shelf Seas in the CMEMS catalogue, http://marine.copernicus.eu/, last access: 27 July 2020).

Figure 1 Environmental data used for the forcing of the dynamic energy budget model (DEB) in the North Sea simulation, showing (a) temperature and (b) chlorophyll a concentration against in situ data from the ICES database.

In the northern Ionian Sea, daily satellite SST data were also obtained from the CMEMS database for the Mediterranean Sea with 0.04^∘ spatial resolution, while CHL a daily data were derived from the Globcolour multi-sensor (i.e. SeaWiFS, MERIS, MODIS, VIIRS, and OLCI a) merged product (http://globcolour.info  last access: 27 July 2020) at 1 km spatial resolution based on the OC5 algorithm suitable for coastal regions (Gohin et al., 2002). The forcing data were averaged over the study area (39.49–39.65^∘ N, 20.09–20.23^∘ E) covering the period 2014–2015 (2 years), when the cultured mussel is ready for the market. The chosen CHL a dataset was found preferable, as compared with other available remote sensing datasets (i.e. CMEMS chlorophyll product for Mediterranean Sea), since it presented a better spatial and temporal coverage (El Hourany et al., 2019; Garnesson et al., 2019) and a slightly lower error, as compared with the very few available in situ data in the study area (not shown). Unfortunately, these were very scarce, and therefore an extended comparison between remote and in situ data could not be conducted. Satellite data have facilitated large-scale ecological studies by providing maps of phytoplankton functional types and sea surface temperature (Raitsos et al., 2005, 2008, 2012, 2014; Palacz et al., 2013; Di Cicco et al., 2017; Brewin et al., 2017). The daily environmental forcing data are shown in Figs. 1 and 2 for the North Sea and the N Ionian Sea respectively. The two coastal environments present some important differences regarding both CHL a and SST. Specifically, in the N Ionian Sea, CHL a is relatively low (annual mean ∼0.88 mg m⁻³) and peaks during winter (maximum ∼2.64 mg m⁻³ at December 2014), while in the North Sea CHL a is about 4 times higher (annual mean 4.25 mg m⁻³), peaking in April every year (maximum range 29.44–33.38 mg m⁻³), as soon as light availability reaches a critical level (Van Beusekom et al., 2009). The higher productivity during the spring season in the North Sea is related with the nutrient inputs from the English Channel, the North Atlantic, and particularly the river discharge of nutrient-rich waters along the Belgian–French–Dutch coastline, which peaks earlier, during winter period (Van Beusekom et al., 2009). The SST peaks during August in both areas (Figs. 1 and 2) but is significantly higher in the N Ionian Sea (maximum 28.8 ^∘C), as compared to the North Sea (maximum 18–19.3 ^∘C).

Figure 2 Environmental data used for the forcing of the dynamic energy budget model in the northern Ionian Sea simulation, showing (a) temperature and (b) chlorophyll a concentration.

The environmental concentration of MPs, C_env (particles L⁻¹), was also obtained at a daily time step as randomly generated values of the Gaussian distribution that is determined by the mean value and standard deviation of the observed field data (0.4±0.3 particles L⁻¹, North Sea; Van Cauwenberghe et al., 2015; 0.0012±0.024 particles L⁻¹, N Ionian Sea, Digka et al., 2018a). Considering that these values originate from surface waters and that mussels live in the near-surface layer (0–5 m), C_env is estimated as a mean value of the upper layer with the methods described by Kooi et al. (2016), which studied the vertical distribution of MPs, considering an exponential decrease with depth. Specifically, in the N Ionian Sea, mussels were collected from a depth up to 3 m (Digka et al., 2018a), while in the North Sea (Van Cauwenberghe et al., 2015), there is no information, and thus a maximum depth of 5 m is adopted. In the North Sea simulation, the effect of tides is taken into account by considering that the mussel originated from the intertidal zone is submerged 12 h during the day (Van Cauwenberghe et al., 2015). In the N Ionian Sea simulation, tides are not considered, given the very small tide amplitude (few centimetres) in the Mediterranean (i.e. Sarà et al., 2011; Hatzonikolakis et al., 2017), and thus the cultured mussel is assumed permanently submerged. In situ hourly tide data (2007–2011) from the coastal zone of the region (Dunkerque station 51.04820^∘ N, 2.36650^∘ E) obtained from Coriolis and Copernicus data provider (http://marine.copernicus.eu, http://www.coriolis.eu.org, last access: 28 July 2020) showed that mussels experience alternating periods of aerial exposure and submergence at approximately every 6 h (two high and two low tides). During aerial exposure, the model suspends the feeding processes (Sarà et al., 2011) and simulates metabolic depression (Monaco and McQuaid, 2018) where the Arrhenius thermal sensitivity equation (Eq. 9) is corrected by a metabolic depression constant (M_d=0.15), a value representative for M. galloprovincialis and here applied also for M. edulis. In the present study, the mussel's body temperature change during low tide is ignored, inducing a model error. The mussel's body temperature (i.e. surrounding water temperature for submerged mussels) during air exposure depends on many factors, such as solar radiation, air temperature, wind speed and wave height, according to studies investigating the temperature effect on intertidal mussels (Kearney et al., 2010; Sarà et al., 2011). However, the present study aims to primarily examine the MPs accumulation, and thus the intertidal mussel's body temperature was not thoroughly examined. Nonetheless, the time that the mussel is able to filter, ingest, and excrete the suspended matter (i.e. food and MP particles) and the effect on the mussel's growth through the modified relation of k(T) are included, since the assimilation process occurs whether the mussel is submerged or not (Kearney et al., 2010).

2.5 Parameter values

Most of the DEB model parameters were obtained from Van der Veer et al. (2006) and refer to the blue mussel M. edulis in the northeast Atlantic (see Table 3 for the exceptions). This assumption has also been adopted in previous studies which showed that this parameter set for M. edulis applies also for M. galloprovincialis (i.e. Casas and Bacher, 2006; Hatzonikolakis et al., 2017). The half-saturation coefficient X_k represents the density of food at which the food uptake rate reaches half of its maximum value and should be treated as a site-specific parameter (Troost et al., 2010; Pouvreau et al., 2006). In order to estimate the value of X_k, a different approach was followed for each study area.

For the North Sea simulation, X_k was tuned so that the simulated individual has the recorded size at the corresponding estimated age (Van Cauwenberghe et al., 2015) growing with the representative growth rates of wild M. edulis at the region (Saraiva et al., 2012; Sukhotin et al., 2007). For the N Ionian Sea simulation, an alternative method was adopted, aiming to generalize the DEB model to overcome the problem of site-specific parameterization. The DEB model was tuned against literature field data for cultured mussels originated from different areas in the Mediterranean and Black seas, where the average CHL a concentration ranged between 1.0 and 5.0 mg m⁻³, and one X_k value was found for each area. The four areas used, their characteristics, and the corresponding value of X_k adopted are shown in Table 4. These values of X_k are related to the prevailing CHL a concentration of each area ([CHL a]) through three different functions: linear: $f\left(x\right)=a\cdot \left[\mathrm{CHL}a\right]+b$; exponential: $f\left(x\right)=a\cdot \exp (b\cdot \left[\mathrm{CHL}a\right])$; and power: $f\left(x\right)=a\cdot {\left[\mathrm{CHL}a\right]}^b+c$. The curve fitting app of MATLAB (MATLAB R2015a) was used for the determination of a,b, and c of each function, taking into account the 95 % confidence level. The score of each function regarding the somatic/mussel growth simulation in all four regions is tested through target diagrams (Jolliff et al., 2009) by computing the bias and unbiased root-mean-square deviation (RMSD) between field and simulated data of all four regions, and the function with the best score is adopted. A similar approach was followed by Alunno-Bruscia et al. (2011) for the oyster Crassostrea gigas in six Atlantic ecosystems, expressing the X_k as a linear function of food density (e.g. phytoplankton). Unfortunately, the approach described for the N Ionian Sea simulation could not be applied in the North Sea, as the limited amount of growth data from the literature for wild M. edulis in similar environments did not permit a statistically significant fit of a similar function, X_k=f(chl a).

Table 4 Half-saturation tuned values (X_k) and mussel growth data (length) in different areas of the Mediterranean and Black seas.

2.6 Simulation of reproduction–initialization of the model

The reproductive buffer (R) is assumed to be completely emptied at spawning (R=0) (Sprung, 1983; Van Haren et al., 1994). In order to simulate mussel spawning, the gonado-somatic index (GSI) defined as gonad dry mass over total dry flesh mass was computed at every model's time step (Eq. 17 Table 1; the water content of the fresh tissue mass was assumed 80 % according to Thomas et al., 2011). Spawning was induced by a critical value of GSI (GSI_th, Table 3) and a minimum temperature threshold (T_th) at each study area, obtained from the literature. In the North Sea implementation, T_th was set at 9.6  ^∘C (Saraiva et al., 2012), while in the N Ionian Sea, at 15 ^∘C (Honkoop and Van der Meer, 1998). This kind of formulation for the spawning event in bivalves has been used in previous studies (i.e. Pouvreau et al., 2006; Troost et al., 2010; Thomas et al., 2011; Monaco and McQuaid, 2018). The simulated abrupt losses of the mussel's tissue mass correspond to spawning events, and the model's prediction was compared with the available literature data regarding the spawning period in each study area. Theodorou et al. (2011) demonstrated that the spawning events occur during winter for M. galloprovincialis in the mussel farms of Greece, while in the North Sea the spawning period for M. edulis is extended from the end of April until the end of June (Sprung, 1983; Cardoso et al., 2007).

In both areas, the model was initialized so that the simulated individual is in the juvenile phase (V<V_p; Table 3) and the reproductive buffer can be considered to be empty (R=0) (Thomas et al., 2011). As stated by Jacobs et al. (2015) amongst others, juvenile mussels (M. edulis) range between 1.5 and 25 mm in size. Specifically, in the North Sea the settlement of mussel larvae (M. edulis) takes place in June and the juveniles grow to a maximum size of 25 mm within 4 months (Jacobs et al., 2014). In the N Ionian Sea, the operating mussel farms follow the life cycle of M. galloprovincialis, starting the operational cycle each year by dropping seed collectors from late November until March, and the juvenile mussels grow up to 6–6.5 cm after approximately 1 year according to the information obtained from the local farms in the region and Theodorou et al. (2011). The initial fresh tissue mass was distributed between the structural volume (V) and reserves energy (E). Energy allocated to those two compartments was firstly constrained by the initial length (L), and then energy allocated to V was in Eq. (10) (Table 1). The initial value of E was set so that the simulated individual has an initial weight that corresponds to the juvenile phase (V<V_p) (Table 5). Finally, for both model implementations, the initial accumulation of MPs in the mussel's tissue (C) was set to zero.

Table 5 Dynamic energy budget-accumulation model: initial values. L: shell length; W: fresh tissue mass; V: structural volume; E: energy reserves; R: energy allocated to reproduction; C: microplastics accumulation.

2.7 Simulation runs

The DEB-accumulation model simulates at an hourly basis the growth and MPs accumulation of the wild mussel from the North Sea and the cultured mussel from the N Ionian Sea. Initially, a model run is performed at each study area during the periods from July 2007 to August 2011 (4 years) for the North Sea simulation and late November 2014 to January 2016 (∼1 year) for the N Ionian Sea simulation. Additionally, the inverse simulations were performed in order to evaluate the depuration phase of both cultured and wild mussel, by setting the environmental MPs concentration equal to zero (C_env=0), after a period of 1 year simulation at the N Ionian Sea, when the cultured mussel has the appropriate size for market, and after 4 years in the North Sea, when literature field data are available (Van Cauwenberghe et al., 2015). In this simulation, the mussel's gut clearance is achieved by the excretion of MPs through faeces (third term of Eq. 18), and thus it is necessary to maintain the existence of food in the mussel's environment in order to ensure that the feeding–excretion processes will occur.

Furthermore, to examine the model's uncertainty related to the environmental MPs concentration, a series of 15 and 13 simulations were performed in the North Sea and N Ionian Sea respectively, adopting different constant values of C_env within the observed range of each area. Finally, the effect of the environmental forcing data and some model parameters on the resulting MPs accumulation by both mussels was explored through sensitivity experiments. These were used to derive a new function that predicts the level of MPs pollution in the environment.

2.8 Sensitivity tests and regression analysis

The effect of the environmental data (CHL a, temperature, C_env) and two parameters representative of mussel's growth (X_k, Y_k) on the MPs accumulation by the mussel for each study area was examined through sensitivity experiments with the DEB-accumulation model. Each variable (CHL a, T, C_env) and parameter (X_k, Y_k) was perturbed by ±10 % to examine its effect on the simulated MPs accumulation, and the results of each run were analysed using a sensitivity index (SI). SI calculates the percentage change in the mussel's MPs accumulation; $\mathrm{SI}=\frac{\mathrm{1}}{n}\sum _{t=\mathrm{1}}^n\frac{\left|C_t^{\mathrm{1}}-C_t^{\mathrm{0}}\right|}{C_t^{\mathrm{0}}}\cdot \mathrm{100}$ (%), where n is the simulated time steps, $C_t^{\mathrm{0}}$ is the MPs accumulation predicted with the standard simulation at time t, and $C_t^{\mathrm{1}}$ is the MPs accumulation with a perturbed variable/parameter at time t; for details see Bacher and Gangnery (2006). The same method has been also applied to other studies, which examined the model's sensitivity for specific variables/parameters regarding the mussel growth (Casas and Bacher, 2006; Rosland et al., 2009; Béjaoui-Omri et al., 2014; Hatzonikolakis et al., 2017). In order to also examine the effect of tides, in the North Sea implementation, the sensitivity experiments were conducted twice as follows: the first time assuming that the mussel is permanently submerged and the second time assuming that the mussel is periodically exposed to the air.

Preliminary sensitivity experiments showed that the MPs accumulation is highly dependent on the prevailing conditions regarding the CHL a, temperature, and C_env and the mussel's growth that is regulated by the half-saturation coefficient (X_k). Therefore an attempt was made using the model's output to describe the MPs accumulation as a function of these variables through a custom regression model as follows:

$$\begin{array}{}\text{(19)}& y=b_{\mathrm{1}}\cdot W+b_{\mathrm{2}}\cdot \exp \left({\displaystyle \frac{\mathrm{1}}{T}}\right)+b_{\mathrm{3}}\cdot {\displaystyle \frac{\mathrm{1}}{\left[\mathrm{CHL}a\right]}}+b_{\mathrm{4}}\cdot C_{\mathrm{env}},\end{array}$$

where y (particles per individual) is the response variable and represent the predicted MPs accumulation by the mussel; W (g) is the mussel's fresh tissue mass, T (K) is the sea surface temperature; CHL a and C_env are the concentrations of chlorophyll a and MPs in the water respectively, which are the predictor variables. The values of coefficients b₁, b₂, b₃, and b₄ are calculated using the nonlinear regression function (nlinfit, MATLAB R2015a), which attempts to find values of the parameters b that minimize the least-squared differences between the model's MPs accumulation output C and the predictions of the regression model $y=f(W,T,$ [CHL a], C_env, b).

The ultimate aim of this analysis, once coefficients are determined, is to use Eq. (19) to obtain the environmental MPs concentration:

$$\begin{array}{}\text{(20)}& C_{\mathrm{env}}={\displaystyle \frac{\mathrm{1}}{b_{\mathrm{4}}}}\cdot \left(C-b_{\mathrm{1}}\cdot W-b_{\mathrm{2}}\cdot \exp \left({\displaystyle \frac{\mathrm{1}}{T}}\right)-b_{\mathrm{3}}\cdot {\displaystyle \frac{\mathrm{1}}{\left[\mathrm{CHL}a\right]}}\right),\end{array}$$

which could be a very useful tool to predict the MPs concentration in the environment, when all involved variables are known (mussel's accumulated MPs, C; wet weight, W; temperature, T; and CHL a), using the mussel as a potential bioindicator (Li et al., 2016, 2019). The score of this custom model was tested by applying Eq. (20) in our study areas and six more areas around the UK, where information on a mussel's wet weight and both the mussels' and environment's MPs load is available (Li et al., 2018). CHL a and temperature, which were not included in Li et al. (2018), were obtained from daily satellite images (same source as in the North Sea; see Sect. 2.4), covering the period that the mussels were harvested (Li et al., 2018).

3 Results

3.1 Growth simulations

The growth simulations of M. edulis and M. galloprovincialis for the North Sea and the N Ionian Sea are shown in Figs. 3 and 4 respectively. In the North Sea implementation, X_k was tuned to a constant value: X_k=8 mg m⁻³ (±1.5 mg m⁻³). The fitted value was higher, as compared to the one (X_k=3.88 mg m⁻³) used by Casas and Bacher (2006) in productive areas of the French Mediterranean shoreline (average CHL a concentration: 1.45 mg m⁻³; maximum peak at 20 mg m⁻³), as a consequence of the higher productivity in the North Sea (average CHL a concentration: 4.25 mg m⁻³; maximum peak at ∼33.40 mg m⁻³). The high value of X_k could also be explained by the presence of inedible particles (i.e. MPs) that led to lower-quality food in the mussel's diet compared with an assumed clean-from-inedible-particles environment (Kooijman, 2006; Ren, 2009). In the present study the inedible particles (i.e. MPs) have been incorporated in the mussel's diet through the modified relation of the functional response f (Eq. 5, Table 1), which regulates the assimilation rate and thus the mussel's growth. However, the DEB model applied at the French site, did not account for inedible particles in the mussel's food. Furthermore, it has been reported that wild mussels grow considerably slower than farmed mussels (∼1.7 times) (Sukhotin and Kulakowski, 1992), and thus, a higher value of X_k promotes less mussel growth, which is the case for the North Sea mussel. The simulated mussel shell length after 4 years, in August, is 4.35 cm, and the fresh tissue mass is 1.87 g, in agreement with Van Cauwenberghe et al. (2015) and other studies conducted on wild mussels (Sukhotin et al., 2007; Saraiva et al., 2012; MarLIN, 2016). In particular, Saraiva et al. (2012) found that after 16 years of simulation, the wild mussel of the Wadden Sea (North Sea) is 7 cm long, while according to Bayne and Worral (1980) a mussel with shell length 4 cm corresponds to the age of 4 years, in agreement with the current study. The simulated growth presents a strong seasonal pattern, being higher during the spring and summer season, as compared to autumn and winter, which is consistent with the seasonal cycle of temperature and CHL a concentration, for a typical year in the region (Fig. 1). The increase in food availability and temperature during spring (April) results in high mussel growth for a 4-month period, while the decrease in CHL a from summer until the end of the year, in conjunction with the temperature decrease in autumn, results in a less mussel growth. Spawning events that occurred in late April–early May (30 April–2 May) each year are responsible for the sharp decline in a mussel's fresh tissue mass, shown in Fig. 4 (Handa et al., 2011; Zaldivar, 2008), which is in agreement with the literature (Sprung, 1983; Cardoso et al., 2007; Saraiva et al., 2012). The predicted weight loss due to spawning was around 7 % at the first year of simulation, while the second, third, and fourth year the percentage of weight loss increased gradually to 8.3 %, 12.6 %, and 14.4 % respectively. Bayne and Worral (1980) demonstrated that the weight losses on spawning for individuals of 1 g weight vary between 2.1 % and 39.8 %, presenting a weight-specific increase with size.

Figure 3 (a) Simulated mussel shell length (L) and (b) fresh tissue mass (W) against North Sea data (red star: mean ± SD), using chlorophyll a, X=[CHL a], in the mussel diet.

Figure 4 (a) Simulated mussel shell length (L) and (b) fresh tissue mass (W) against North Sea data (red star: mean ± SD), using chlorophyll a, X=[CHL a], in the mussel diet.

In the N Ionian Sea implementation, X_k is applied as a function of CHL a concentration through the method described in Section 2.5. The target diagram showing the performance of each tested function (linear: $f\left(x\right)=a\cdot \left[\mathrm{CHL}a\right]+b$, where a=0.959 and $b=-\mathrm{1.420}$; exponential: $f\left(x\right)=a\cdot \exp (b\cdot \left[\mathrm{CHL}a\right])$, where a=0.2 and b=0.567; power: $f\left(x\right)=a\cdot {\left[\mathrm{CHL}a\right]}^b+c$, where a=0.01, b=3.529 and c=0.480) is shown in Fig. 5. The linear and power function of X_k present good skill, with the power function leading to the most successful simulation of the cultured mussel's growth in all four areas (diagram marks for mussel length and fresh tissue mass are closer to the target's centre). The power function applied in the N Ionian Sea resulted in a mussel's shell length of 5.8 cm and fresh tissue mass of 5.92 g after 1 year of simulation, in agreement with Theodorou et al. (2011). The spawning event occurred at the beginning of December (Theodorou et al., 2011) and was illustrated by a 12.6 % tissue mass decline.

Figure 5 Target diagram of simulated shell length (L) and fresh mass tissue weight (W) against field data from the Thermaikos and Maliakos gulfs (eastern Mediterranean Sea), Black Sea, and Bizerte Lagoon (southwestern Mediterranean Sea), using the power (L₁,W₁), exponential (L₂,W₂), and linear (L₃,W₃) functions of the half-saturation coefficient. The model bias is indicated on the y axis, while the unbiased root-mean-square deviation (RMSD) is indicated on the x axis.

3.2 Microplastics accumulation and depuration phase

The hourly simulated MPs accumulation by the mussel in the North Sea and N Ionian Sea are shown in Figs. 6 and 7 respectively. Calibration of the parameter k_f (1.2 d⁻¹) led to a model which was well fitted to the observed MPs accumulation in the mussel of both study areas. In the North Sea, a 4-year-old wild mussel (L=4.35 cm, W=1.87 g) contains 0.53 particles per individual in August, within the range value found by Van Cauwenberghe et al. (2015) (0.4±0.3 particles per individual), although the model overestimated the data range, reproducing a seasonal increase that was not observed. This is most likely due to the fact that Van Cauwenberghe et al. (2015) allowed a 24 h clearance period before analysing the mussels' tissue for MPs, resulting in slightly lower MPs accumulation than the model's prediction. The MPs egested through faeces by the 4-year-old mussel after 24 h were 0.2±0.2 particles per individual (Van Cauwenberghe et al., 2015), which agree also with model's output (0.3 particles per individual, Fig. 8) regarding the depuration phase and could compensate for the observed difference in the mussel's MPs load between the simulated and field data. In the N Ionian Sea, the simulated MPs accumulation by the cultured mussel with L=4.85 cm and W=3.33 g was 0.91 particles per individual at the end of June, in agreement with field observations obtained from Digka et al. (2018a) (0.9±0.2 particles per individual). Overall, the developed model simulated the MPs accumulation by both mussels in the two different areas, using the same parameter set (see Table 3 for the exceptions), under the assumption that parameters referred to as silt particles (i.e. inedible particles) may also be used to describe the MPs accumulation. Both simulations were in good agreement with the available field data, with a small deviation for the North Sea. This may lead to the assumption that mussels present a common behaviour against all inedible particles. In the model's results, based on the uptake and excretion rates of MPs by the mussels in both study areas, the majority of MPs are rejected through pseudofaeces and fewer through faeces production (not shown). This is in agreement with Woods et al. (2018), which found that most microplastic fibres (71 %) were quickly rejected as pseudofaeces and <1 % excreted in faeces.

Figure 6 Microplastics (MPs) accumulation by the mussel (blue line) against field data (red star: mean ± SD), using daily environmental concentration of MPs (C_env mean value ± SD: 0.4±0.3 particles L⁻¹) in the North Sea.

Figure 7 Microplastics (MPs) accumulation by the mussel (blue line) against field data (red star: mean value ± SD), using daily environmental concentration of MPs (C_env mean value ± SD: 0.0012±0.024 particles L⁻¹) in the northern Ionian Sea.

The small-scale (daily) fluctuations of MPs in the mussel (wild and cultivated) reflect the adopted random variability in the environmental MPs concentration C_env and the daily fluctuations of the environmental forcing (CHL a, temperature). The large-scale (seasonal) variability follows mainly the variability of the clearance rate. The seasonal variability in the CHL a concentration and temperature greatly determines the variability in the clearance rate and hence the variability in MPs in the individual. Moreover, the model predicts that mussel's energy needs are increased as it grows, and therefore the clearance rate is increased, resulting in higher MPs accumulation.

The simulated time needed to clean the mussel's gut from the MPs load for both areas is shown in Fig. 8. In both areas, the cleaning follows an exponential decay, in agreement with laboratory experiments by Woods et al. (2018). In particular, the model predicts a 90 % mussel's cleaning after 284 h (∼12 d) and 56 h (∼2.5 d) for the N Ionian Sea and North Sea respectively. The cleaning process is more rapid in the North Sea simulation, which can be attributed to the higher CHL a concentration found in this area, leading to increased production of faeces by the mussel and hence faster excretion of the accumulated MPs. In the N Ionian Sea, on the other hand, the rate of the mussel's cleaning is slower, due to the limited food availability.

Figure 8 Depuration phase of the cultured Mytilus galloprovincialis (red line) and wild Mytilus edulis (black line) using zero environmental concentration of microplastics (C_env=0) after 1 and 4 years of simulation time in the northern Ionian Sea and North Sea respectively.

3.3 Model's uncertainty regarding the environmental microplastics concentration

The MPs concentration in the environment presents a strong variability in both temporal and spatial scales. To examine the model's uncertainty related to the environmental MPs concentration (C_env), a series of 15 and 13 simulations were performed in the North Sea and N Ionian Sea respectively, adopting different values of C_env within the observed range of each area. In the North Sea, the adopted C_env ranged between 0.1 and 0.8 particles L⁻¹ with a step of 0.05 (15 runs), while in the N Ionian Sea C_env ranged between 0.0012 and 0.0252 particles L⁻¹ with a step of 0.002 (13 runs). The mean seasonal values and standard deviation of the 15 simulations in the North Sea and the mean monthly values and standard deviation of the 13 simulations in the N Ionian Sea were computed and plotted in Figs. 9 and 10 respectively. Each error bar represents the uncertainty in the simulated accumulation at the specific time, related to the environmental MPs concentration.

Figure 9 Mean seasonally values and standard deviation of microplastics (MPs) accumulation (red error bars: mean value ± SD) by the mussel in North Sea derived from 15 model runs with different constant values of environmental MPs concentration (C_env range: 0.1–0.8 particles L⁻¹); Mean hourly simulated data (black line) and standard deviation (blue lines) of microplastics accumulation derived from three model runs with stochastic sequences of daily random C_env values.

Figure 10 Mean monthly values and standard deviation of microplastics accumulation (red error bars: mean value ± SD) by the mussel in northern Ionian Sea derived from 13 model runs with different constant values of environmental MPs concentration (C_env range: 0.0012–0.024 particles L⁻¹); Mean hourly simulated data (back line) and standard deviation (blue lines) of microplastics accumulation derived from three model runs with stochastic sequences of daily random C_env values.

In both case studies, the uncertainty of the model appears to increase as the MPs accumulation is increased. As the mussel grows in the North Sea, the mean value and standard deviation of MPs accumulation is increased during the same season every year, illustrating the effect of the mussel's weight. Moreover, the seasonal variability in the MPs accumulation appears to be related with the seasonality of CHL a concentration. This is apparent during each year's spring, when CHL a concentration peaks at its maximum value (∼30 mg m⁻³; see Fig. 1); the filtration rate is decreased (Riisgård et al., 2003, 2011), leading to lower MPs accumulation by the mussel and thus lower model uncertainty. In the N Ionian Sea, the effect of the mussel's weight is more apparent in the early months (∼6 months), resulting in higher MPs accumulation and model uncertainty as the mussel grows. Afterwards, the seasonality of both CHL a concentration and temperature plays the major role. During summer, when the CHL a concentration is progressively decreased, reaching minimum values (∼0.7 mg m⁻³), and temperature is increased (>20 ^∘C), the filtration rate is significantly decreased or stopped, resulting in lower MPs accumulation and lower model uncertainty. This is in line with studies reporting that the mussel suspends the filtering activity and thus closes its valves until better conditions occur (Pascoe et al., 2009; Riisgård et al., 2011). Overall, the available field data lie within the model's uncertainty, apart from the North Sea case, where the range of field data variability and model uncertainty dot not overlap significantly at the time of the observations.

Moreover, to evaluate the scenario adopted with the set-up of the previous experiments (random C_env at a daily time step) three additional model runs are performed in each study area, adopting each time different stochastic sequences of daily random C_env values within the observed range, which is considered to reflect the high spatial and temporal variability in the environmental MPs concentration. The mean value and standard deviation of these “stochastic” runs lie most of the time within the standard deviation of the overall model uncertainty in both case study areas (Figs. 9 and 10).

3.4 Sensitivity and regression analysis results

The results of the sensitivity experiments regarding the MPs accumulation by the mussels are shown in Figs. 11 and 12 for the North Sea and N Ionian Sea respectively. The comparison between the intertidal and subtidal mussel of the North Sea revealed that both +10 % and −10 % perturbation of CHL a and X_k have a slightly lower effect on the MPs accumulation by the intertidal mussel, which is probably attributed to the intermittent feeding periods experienced by the individual due to the tide effect. As far as the temperature effect, both a +10 % and −10 % perturbed value led to higher sensitivity to the MPs accumulation by the intertidal mussel, due to the adopted modified temperature relation during low tide. Especially, if the mussel's body temperature change during air exposure would be considered, the perturbed temperature will probably affect the MPs accumulation even more for the intertidal than the subtidal mussel. The sensitivity of the C_env to the MPs accumulation when perturbed either +10 % or −10 % is almost the same for the intertidal and subtidal mussel, indicating that the environmental MPs concentration affects similarly both mussels, regardless the continuous or intermittent feeding–excretion process.

Figure 11 Sensitivity index of MPs accumulation for the wild mussel of the North Sea when variables (CHL a, temperature, C_env) and parameters (X_k,Y_k) are perturbed ±10 %. The notation (s) refers to the permanently submerged mussel.

Figure 12 Sensitivity index of MPs accumulation for the cultured mussel of the northern Ionian Sea when variables (CHL a, temperature, C_env) and parameters (X_k,Y_k) are perturbed ±10 %.

The comparison between the mussel sensitivity indexes in the N Ionian and the North Sea (in conditions of submergence) study areas reveals some important differences. Generally, most of the perturbed (either +10 % or −10 %) variables and parameters (i.e. CHL a, temperature, X_k) present higher sensitivity to the MPs accumulation by the mussel from the N Ionian Sea. This is attributed to the prevailing environmental conditions and specifically the lower food availability (CHL a) and the higher temperature range in the N Ionian Sea compared to the North Sea, which greatly determine the feeding processes, the mussel's growth, and hence the MPs accumulation. The perturbed C_env in both study areas appears to affect similarly the MPs accumulation for both mussels (∼10 %), with the small difference (<2 %) probably attributed to the higher abundance of seawater MPs present in the North Sea compared to the N Ionian Sea. Finally, the half-saturation coefficient for the inorganic particles (Y_k) has no effect on the MPs accumulation of both North Sea and N Ionian Sea mussels, indicating that the amount of inedible particles (i.e. MPs) is relatively low in both areas, and thus the Y_k does not affect the way that the organic particles are being ingested (Kooijman, 2006). According to Ren (2009), when the inorganic matter is low, the K(y) (Eq. 5; Table 1) is approximately equal to X_k, and then Y_k is the least sensitive parameter for the ingestion rate and thus growth.

The DEB-accumulation model output was used to determine the coefficients in Eq. (19) by the nonlinear regression analysis: b₁=0.1909 (±0.0006), b₂=0.0412 (±0.0019), b₃=0.1315 (±0.0021), and b₄=1.1060 (±0.0253). The accurate estimation of the coefficients (b₁, b₂, b₃, b₄) is indicated by the low confidence intervals, while the mean squared error of the regression model appears also sufficiently small (MSE =0.0523). Subsequently, as shown in Fig. 13, Eq. (20) may be used to predict the MPs concentration of the environment where mussels live. In most cases, the predicted MPs concentration is found within the standard deviation of the field data. Two exceptions are shown in Hastings-A and Plymouth areas. The reasons behind these discrepancies may be related to the environmental conditions prevailing in each area at the sampling time. For example, Eq. (20) does not take into account the impact of tides that may affected the mussel's MPs load (C), and the lack of information on the exact sampling date led to using a mean SST and CHL a value representative of the given sampling time period (Li et al., 2018). Although Eq. (20) does not account for the tide effect, the sensitivity analysis (Fig. 11) showed that the effect of C_env on the mussel's MPs accumulation was the same for both the intertidal and subtidal mussel in the North Sea. This result may also apply at the two exceptions areas, leading to the assumption that the discrepancies are due to the lack of the ambient temperature and CHL a information during the sampling date. In any case, this is a first rough demonstration of the method and should be implemented in more environments in order to be further validated.

Figure 13 Prediction of seawater microplastics concentration by applying Eq. (20) in the northern Ionian Sea, North Sea (present study), and six areas around the UK (Filey, Hastings-A, Hastings-B, Brighton, Plymouth, and Wallasey; Li et al., 2018).

4 Discussion

A DEB-accumulation model was developed and validated with data available from the North Sea and the N Ionian Sea, to study the MPs accumulation by wild M. edulis and cultured M. galloprovincialis, grown in different, representative environments. Although the study is limited by scarce validation data, it should be noted the MPs accumulation model parameter set, except one tuning parameter (k_f), was extracted from the literature (Table 3), assuming that mussels adopt a common defensive mechanism against inedible particles (i.e. silt, MPs). Thus, the theoretical background constructed by Saraiva et al. (2011a) (based on Kooijman, 2010) regarding the feeding and excretion processes of the mussel remains unspoiled. Through the strong theoretical background of DEB theory, this study highlights that the accumulation of MPs by the mussel is highly dependent on the prevailing environmental conditions which control the amount of MPs that the mussel filtrates and excretes.

Towards a generic DEB model, the applied function of the half-saturation coefficient, $f\left(x\right)=a\cdot {\left[\mathrm{CHL}a\right]}^b+c$, successfully captures the physiological responses and thus the growth rate of the cultured mussel at the N Ionian Sea implementation. In the current study, this method led to a robust and generic DEB growth model able to simulate the mussel growth in representative mussel habitats of the Mediterranean Sea, covering a range of productivity and sea surface temperature. This approach supports and takes one step further by using the Bourlès et al. (2008) suggestion about a seasonally varied half-saturation coefficient, demonstrating an improvement of the food quantifier. The applied function of X_kconsiders the daily CHL a fluctuations and, thus, the seasonal variation of the seawater composition. As more field data become available from various environments, the applied approach could result to more generic formulations for the site-specific parameter X_k, so that the model could be applied in several areas of interest, where field growth data are absent and/or to simulate the potential mussel growth in the 2D space.

The simulation of MPs accumulation by the mussels, using the DEB-accumulation model, is in good agreement with the available field data (Figs. 3 and 4). The simulated values lie within the observed field data range (mean ± SD), although the seasonal increase reproduced by the model in the North Sea implementation did not exactly overlap with the field data at the time of observations. This could be attributed to the clearance period (24 h) that allowed mussels to excrete MPs through faeces (0.2±0.2 particles per individual) before the mussel's tissue analysis (Van Cauwenberghe et al., 2015). The measured loss of the mussel's MPs is in agreement with the model's result for the depuration experiment after 24 h. The MPs accumulation by the cultivated mussel (fresh tissue mass 3.33 g) originated from the N Ionian Sea with mean $C_{\mathrm{env}}=\mathrm{0.0012}\pm \mathrm{0.024}$ particles L⁻¹ is 0.91 particles per individual and by the wild mussel (fresh tissue mass 1.87 g) from the North Sea with mean $C_{\mathrm{env}}=\mathrm{0.4}\pm \mathrm{0.3}$ particles L⁻¹ is 0.53 particles per individual. If these concentrations are expressed per gram of wet tissue of mussel, the cultivated mussel contamination (0.27 particles g⁻¹ w.w.) is comparable with the wild mussel (0.28 particles g⁻¹ w.w.), despite the much lower environmental MPs concentration (C_env) in the N Ionian Sea than in the North Sea. This comparison aims to highlight the significant impact of the prevailing environmental conditions (CHL a and temperature) on the MPs accumulation by the mussels, although they originate from different areas and lived during different time periods. The generally high abundance of CHL a in the North Sea simulation contributes to a reduction of the filtering activity and hence of the MPs accumulation. The threshold algal concentration for reduction of the mussel's filtration rate (incipient saturation) has been found to lie between 6.3 and 10.0 mg m⁻³ (Riisgård et al., 2011), which is a range comparable to the CHL a concentrations in the North Sea. Van Cauwenberghe and Janssen (2014) found that cultivated M. edulis from the North Sea contained on average 0.36±0.07 particles g⁻¹ w.w., a slightly higher value than that found in the present study for the wild mussel of the North Sea (0.28 particles g⁻¹ w.w.). This could be attributed to mussel farms acting as a potential source of MPs contamination for the mussels due to plastic materials (i.e. plastic sock nets and polypropylene long lines) used during cultivation (Mathalon and Hill, 2014; Santana et al., 2018). Moreover, the intertidal wild mussel (present study) is assumed to filter and excrete MPs half of the time in comparison with the submerged cultured mussel in the North Sea, resulting though in similar accumulation level. The model also predicts the time needed for the 90 % gut clearance of both cultured (N Ionian Sea) and wild (North Sea) mussels to be almost 284 and 56 h (equivalent to 12 and 2.5 d) respectively, when MPs contamination is removed from their habitat. This is in line with a series of studies which demonstrated that the depuration time varies between 6–72 h and can last up to 40 d depending on several factors such as species, environmental conditions (Bayne et al., 1987), size, and type of MPs (Browne et al., 2008; Ward and Kach, 2009; Woods et al., 2018; Birnstiel et al., 2019).

The strong dependence of food (CHL a), temperature, and seawater MPs concentration on the MPs accumulation by the mussel, regarding its wet weight, is demonstrated through sensitivity experiments that were used to derive a rather simple nonlinear regression model (Eq. 19). The comparison of the regression model's with the DEB model's output resulted in a quite accurate estimation of the coefficients, which in turn sparked the idea of a “new” relationship (Eq. 20) that could potentially predict the MPs concentration in the environment (C_env) when certain conditions are known (CHL a, T, C, W). The latter equation was applied in eight areas in total (two from the present study areas and six from Li et al., 2018), with relatively good results since there is general overlapping of regressed and observed MPs concentration in the environment (C_env), except for Hastings-A and Plymouth areas, probably due to missing information on the environmental conditions (CHL a, SST) during the sampling, suggesting that the mussels can be used as potential bioindicators. Mussels have been previously proposed as bioindicators for marine microplastic pollution (<1 mm), although the efficient gut clearance and selective feeding behaviour limit their quantitative ability (Lusher et al., 2017; Bråte et al., 2018; Beyer et al., 2017; Fossi et al., 2018; Li et al., 2019).The recent study by Ward et al. (2019b) demonstrated that bivalves are poor bioindicators of MPs pollution due to the particle selection during feeding and excretion processes that is based on the physical characteristics of the MPs. Considering that the MPs accumulation is site-dependent and that sampling of mussels is usually easier than seawater (Karlsson et al., 2017; Bråte et al., 2018), models like the one described in Eq. (20), besides the MPs accumulation, take into account also characteristics of the environment that are crucial for the way that mussels accumulate MPs. This method could be possibly used at a global level and allow comparisons between various environments. However, the method described should be validated in more environments with more frequent field data to be able to provide secure results.

In addition to the scarce validation data regarding the MPs accumulation in mussels, this study has some more limitations. First of all, the data regarding the concentration of MPs in the mussels' environment are also scarce; since MPs are a relatively recent subject of study, the existing knowledge of the spatial and temporal distribution is still quite limited (Law and Thompson, 2014; Browne, 2015; Anderson et al., 2016; de Sá et al., 2018; Smith et al., 2018; Troost et al., 2018). To overcome the lack of environmental MPs time series, a function of randomly generated values within the observed range of each area was applied, and its uncertainty was examined through an ensemble forecasting. Specifically, the model's uncertainty due to the environmental MPs concentration (C_env) was tested by performing a series of model runs forced by an envelope of representative values of C_env. The results (Sect. 3.3) showed that the adopted stochastic scenario simulated quite satisfactorily the MPs accumulation by the mussels, lying within the observed field range, although a slight overestimation was found in the North Sea. The approach used is assumed to represent the natural variability since it has been reported that tides, wind, wave action, ocean currents, river inputs, and hydrodynamic features lead to high spatial and temporal variability in MPs distribution even on very small scales (Messinetti et al., 2018; Goldstein et al., 2013). In addition, the nature of the variable C_env makes it difficult to estimate, presenting large observational errors, not only due to the intense physical variation but also due to different sampling and analysis techniques that were used. In a future work the DEB-accumulation model could be coupled with a high-resolution MPs distribution model (Kalaroni et al., 2019), being extensively validated against field data that will have been collected and processed according to a common scientifically defined protocol, to overcome this limitation. Moreover, the approach followed in calculating the value of MPs concentration in the near-surface layer (0–5 m depth) (Kooi et al., 2016) resulted in a representative value of the upper-ocean layer. In-depth knowledge of the MPs distribution, both horizontally and vertically, is essential to understand and mitigate their impact not only on the various marine compartments but also on the organisms inhabiting those compartments (Van Sebille et al., 2015; Kooi et al., 2016). For that reason, it is important to enhance the monitoring activity especially in the vulnerable coastal environments, adopting integrated cross-disciplinary approaches and monitoring of biological, physical, and chemical parameters which provide information on the ecosystem function, in order to improve the assessment of emerging pollutants (i.e. MPs) and their impacts on biota (objective of JERICO-RI framework).

Our assumption that the mussel has the same filtration rate for all particles, independent of their chemical composition, size, and shape, is a simplification and an open theme of discussion (see Saraiva et al., 2011a for details). However, in our model application, a pre-ingestive particle selection by the mussel is implied based on the organic–inorganic content of the suspended matter illustrating the different binding probabilities applied for algal and MP particles during the ingestion process. Through an investigation of wild mussel's faeces and pseudofaeces production in laboratory conditions, Zhao et al. (2018) found that the length of MPs was significantly longer in pseudofaeces than in the digestive gland and faeces. Furthermore, Van Cauwenberghe et al. (2015) demonstrated that mussel's faeces contained larger MPs (15–500 µm) compared to the mussel's tissue (20–90 µm). Apparently, smaller-sized MPs seem to be dominant within the mussels in comparison with the size of the MPs in the ambient environment (Li et al., 2018; Qu et al., 2018; Digka et al., 2018b), implying that the mussel is more prone to ingest and retain smaller-sized MPs. As an example, Digka et al. (2018b) confirmed that the smaller MPs (<1 mm) occupy 62.3 %, 96.9 %, and 100 % of the total MPs in seawater, sediments, and mussels from the N Ionian Sea respectively. In a future work this selection pattern regarding size could be simulated by suitable preference weights among different MPs sizes. This will improve the knowledge of the feeding and excretion mechanisms used by the mussels against MPs pollution and the assessment of the ecological footprint (Rist et al., 2019).

Our assumption that the contamination by MPs does not affect the energy budget in terms of growth might also be a simplification as this is a subject currently under investigation. Van Cauwenberghe et al. (2015) found that although mussels M. edulis exposed to MPs increased their energy consumption, the energy reserves were not affected compared to the control organisms, implying that mussels are able to adopt a defensive mechanism against the suspended inorganic particles (i.e. MPs) (Ward and Shumway, 2004). Furthermore, MPs exposure showed no significant effect on a mussel's (Perna perna) energy budget, despite its long duration and relatively realistic intensity, leading to the hypothesis that mussels can acclimate to the MPs exposure to maintain their health (Santana et al., 2018). On the contrary, other authors suggested a significant energy shift from reproduction to structural growth and elevated maintenance costs, probably attributed to the reduced energy intake, when the organisms (i.e. oyster Crassostrea gigas) were contaminated with high and unrealistic concentration of MPs (Sussarellu et al., 2016). Moreover, Gardon et al. (2018) showed that the overall energy balance of oyster Pinctada margaritifera was significantly impacted by the reduced assimilation efficiency in correlation with the exposed dose of MPs, and for that reason energy had to be withdrawn from reproduction to compensate for the energy loss. In the future, dedicated experiments exploring the effects on all components of a DEB model should be carried out considering long-term realistic MPs exposure.

Our use of the tide data led to some model bias, since the model does not take into account the mussel's body temperature change when this is exposed to air. Assessing the mussel's body temperature requires extended experiments in field conditions (Tagliarolo and McQuaid, 2015; Monaco and McQuaid, 2018). The study by Seuront et al. (2019) along the French coast of the eastern English Channel found no significant correlation between air and a mussel's body temperature but demonstrated a significant positive correlation between the body temperature and the hard substrate (i.e. rocks) temperature. However, in the present study the tide effect on processes that are affected by the thermal equation, k(T), is considered indirectly through the metabolic depression (details in Sect. 2.4). Sarà et al. (2011) coupled a DEB model with a biophysical model (Kearney et al., 2010), incorporating the change of mussel's body temperature during emersion by using information of various climatological variables (i.e. solar radiation, air temperature, wind speed, wave height), but ignored the temperature sensitivity for the physiological processes. In a future study, a combined approach of coupling the present DEB-accumulation model with a biophysical model, which includes both the tide effect on the physiological processes and the mussel's body temperature respectively, could be followed and lead to a more detailed simulation of the intertidal mussel.

5 Conclusions

In a future study the model should be corroborated further by using a larger dataset of MPs accumulation, with sampling of mussels of various sizes and life stages. Currently, the model is mainly limited by the insufficient validation, as a larger dataset could be also used for a better model calibration. However, this study provides a new approach in studying the accumulation of MPs by filter feeders and reveals the relations between characteristics of the mussel's surrounding environment and the MPs accumulation, which is presented with high seasonal fluctuations. Additionally, in a future study the DEB-accumulation model will be coupled to a hydrodynamic–biochemical model (e.g. Petihakis et al., 2002, 2012; Triantafyllou et al., 2003; Tsiaras et al., 2014; Ciavatta et al., 2019; Kalaroni et al., 2020) and a MPs distribution model (Kalaroni et al., 2019) that will provide fields of temperature, food availability, and MPs concentration respectively at the Mediterranean scale and eventually lead to an integrated representation of the MPs accumulation by mussels (Daewel et al., 2008). This fully coupled model will be downscaled to the Cretan Sea SuperSite (a marine observatory dedicated to multiple in situ observations at appropriate spatio-temporal resolution, in a restricted geographical region, maintained over long timescales, and designed to address interdisciplinary objectives, driven by science and society needs), while the parameterization of important biological processes will be redesigned based on the new data which will be acquired in the framework of the JERICO S3 project (http://www.jerico-ri.eu, last access: 28 July 2020). The present study highlights the urgent need for adopting a multidisciplinary monitoring activity by measuring physical, biological, and chemical parameters that are crucial for mapping the MPs distribution, assessing the contamination level of the marine organisms, and investigating the impact on the health status. Overall, despite the limitations mentioned, taking into account that plastics are one of the global hot issues, this particular study could help design next efforts, since it provides indications of the future related priority issues.