Northern Europe’s suitability for offshore European flat oyster (Ostrea edulis) habitat restoration based on population dynamics

Stechele, Brecht; Barbut, Leo; Lacroix, Geneviève; van Duren, Luca A.; Van Lancker, Vera; Degraer, Steven; Gavazzi, Giacomo Montereale; Bossier, Peter; Declercq, Annelies M.; Nevejan, Nancy

doi:10.3389/fmars.2023.1224346

ORIGINAL RESEARCH article

Front. Mar. Sci., 07 September 2023

Sec. Marine Ecosystem Ecology

Volume 10 - 2023 | https://doi.org/10.3389/fmars.2023.1224346

Northern Europe’s suitability for offshore European flat oyster (Ostrea edulis) habitat restoration based on population dynamics

Leo Barbut²

Luca A. van Duren³

Giacomo Montereale Gavazzi²

Peter Bossier¹

Annelies M. Declercq¹

Nancy Nevejan¹

¹Department of Animal Sciences and Aquatic Ecology, Ghent University, Gent, Belgium
²Operational Directorate Natural Environment, Royal Belgian Institute of Natural Sciences, Brussel, Belgium
³Deltares, Delft, Netherlands

Introduction: European and member state legislation encourage the sustainable development of offshore environments, with the restoration or creation of flat oyster (Ostrea edulis) habitats offering potential solutions to enhance biodiversity and provide ecosystem services. The site selection and identification of suitable habitats for flat oyster habitat restoration projects remains a significant knowledge gap.

Methods: By combining seabed substrate information with a coupled population (Dynamic Energy Budget - Individual-Based Models) and particle tracking model (for larvae dispersal), critical insights can be gained into the spatial distribution of suitability indicators such as population growth, fitness, reproduction, and self-recruitment. By applying this model to the English Channel and the North Sea over a ten-year period, suitable locations for flat oyster habitat restoration, restorative aquaculture, or oyster-related nature-inclusive designs can be identified. Comparing historical oyster bed locations with model outputs offers insights into why these locations were suitable for oyster bed development and allows for validation of the proposed approach.

Results and discussion: Coastal and nearshore environments are generally more suitable for flat oyster habitat restoration, and populations will grow more quickly in these areas. Offshore restoration in the North Sea presents challenges if it relies solely on self-recruitment. Besides site selection purposes, the model can be used to evaluate the effect of management strategies (e.g., initial population size) or environmental pressures (e.g., climate change, pollution) on restorations success.

1 Introduction

European (EU) frameworks, such as the Marine Strategy Framework Directive, the EU Green Deal, the EU Biodiversity Strategy 2030 and the upcoming Nature Restoration Law, mandate EU member states to protect and restore marine biodiversity, habitats and ecosystems. Additionally, nature restoration and protection are becoming key prerequisites for tenders in the offshore energy and marine infrastructure sector, thereby commercially incentivizing nature restoration and creation.

Oyster beds and reefs provide essential ecosystem functions within coastal and offshore ecosystems around the world, providing structural complexity leading to numerous ecosystem services, as well as contributing to coastal economies for centuries (Beck et al., 2011). In Europe, flat oysters (Ostrea edulis) have been identified as a priority species for protection and habitat restoration (Haelters and Kerckhof, 2009) due to their potential to engineer habitats (Todorova et al., 2009; Kregting et al., 2020), and promote biodiversity (Christianen et al., 2018; Lown et al., 2021; Hemraj et al., 2022). As such, they are the focus of several flat oyster nature-based solution projects, restorative flat oyster aquaculture projects (Bertolini and Pastres, 2022; H2020 UNITED, Horizon ULTFARMS) and a multitude of habitat creation or restoration projects in many European countries. (https://noraeurope.eu).

Upscaling of oyster bed restoration efforts has been difficult and large-scale projects have only been on-going in the USA (Schulte et al., 2009), Australia (Gillies et al., 2015), and Europe (Pogoda et al., 2020). These are mainly situated inshore, in coastal or protected areas, although Europe’s largest beds were historically located offshore (Houziaux et al., 2008; Bennema et al., 2020). Implementation of offshore oyster habitat restoration is limited due to high operational costs and logistical and legal issues, such as access to restricted areas (e.g., wind farms). Almost all locations where historical beds were present, including most European MPAs with a designated Natura 2000 habitat directive status, are regularly impacted by bottom trawling and hence unsuitable for oyster restoration efforts. However, offshore renewable energy projects have successfully excluded bottom disturbing activities from the concession zones to protect their underwater infrastructure. Consequently, offshore wind farms have received a status as quasi-marine protected areas (Hammar et al., 2015) and are considered suitable areas for large scale offshore nature restoration projects including flat oyster bed restoration (Kamermans et al., 2018). Nevertheless, how to evaluate the suitability of the offshore environment for flat oyster bed restoration remains a key knowledge gap (zu Ermgassen et al., 2020).

To this end, the Native Oyster Restoration Alliance (NORA) established the Site Selection Working Group to support flat oyster restoration projects. In response to questionnaires, the NORA Site Selection Checklist (Hughes and zu Ermgassen, 2021; Hughes et al., 2023) was produced and aggregated factors relevant for flat oyster site suitability studies into essential and desirable categories, each of which was further divided into five thematic categories including (1) threats to habitat restoration, (2) project logistics, (3) socioeconomic factors, (4) biotic factors and (5) abiotic factors. While most factors can be addressed through management strategies, biotic and abiotic factors constitute fundamental requirements that are, in most cases, impossible to alter. Consequently, understanding these biotic and abiotic factors serves as a prerequisite for the development of niche models that can be utilized in site selection for these restoration projects.

The spatial allocation of oyster habitat restoration projects is commonly achieved through assignment of restoration efforts in the vicinity of local beds, or organization of restoration efforts on sites where beds historically existed. However, historical absence-presence data, such as the Olsen fisheries map (Olsen, 1883) are often inaccurate for offshore habitats (Bennema et al., 2020), and the environment has changed significantly over the last hundreds of years (Bindoff et al., 2019). Furthermore, historical beds often do not coincide with areas designated as Marine Protected Areas (MPAs) or locations for offshore wind which provide adequate protection from human disturbances. When envisioning oyster habitat creation instead of habitat restoration, one cannot solely rely on historical absence presence data, but should include ecological knowledge in the site selection process.

A fundamental niche model (species-environment interactions) for the European flat oyster has been developed by Stechele et al. (2022) and is defined as a set of sediment type, temperature, salinity, and food availability. Fundamental niche models for individuals can evaluate if the environment supports the individual’s life cycle (Kearney et al., 2010; Hale et al., 2019) including settlement, survival, growth, and reproduction (Glover et al., 2021; Stechele et al., 2022; Stechele et al., 2023). Other environmental factors such as oxygen, turbidity, current speed, or wave action also define the fundamental niche but the impact of these variables on the metabolism of flat oysters is insufficiently understood. However, it is likely that several of these factors are key in determining where oysters can proliferate and where not.

When performing site selection for restoration, one should understand the suitability of a location for flat oyster populations, more than only for flat oyster. By integrating self-recruitment data to population models, it becomes clear at which locations populations will thrive and where restoration efforts are can become self-supportive.

This work establishes a realistic niche model (which include species-environment interactions and in addition, inter and intra species interactions) for flat oyster habitat restoration and applies the model to the English Channel and the North Sea to indicate suitable areas for offshore flat oyster habitat restoration.

2 Methods

2.1 Geographical scale and environmental data layers

The research domain encloses the North Sea and the English Channel 48.5° N to 57° N (resolution 0.0417° N) and -4° E to 9° E (resolution 0.0833° E). The area encompasses some important historical offshore flat oyster beds and restoration projects, and partly corresponds to the geographical range of the North Atlantic offshore oyster populations (Launey et al., 2002). Environmental drivers (chlorophyll a, temperature, salinity) were extracted from Copernicus Marine Service (European Union-Copernicus Marine Service, 2020a; European Union-Copernicus Marine Service, 2020b) for the geographical range of interest with a daily time step. Depth profiles of the environmental drivers were converted to bottom values by extraction of the deepest numerical value. Products were daily 25-hour, de-tided, averages from 1/1/2000 to 31/12/2009. The horizontal and vertical resolution of the products (0.111° × 0.067°) were interpolated to the resolution of the research domain used in this study using a nearest-neighbour interpolation algorithm. The Copernicus Marine Environment Monitoring Service (CMEMS) products do not cover inshore or estuarine locations and therefore this study is only applicable to nearshore (0 – 20 km from the shore) and offshore (>20 km from the shore) locations. European Marine Observation and Data Network (EMODnet) multiscale substrate maps (combination of 250k and higher resolutions, where available) (EMODnet Geology, 2023) were extracted. Folk 5 classification, geographically most widely available, was chosen to depict the distribution of mud to muddy sands, sands, mixed sediments, coarse substrate, as well as rock and boulders (Kaskela et al., 2019). Locations of historical beds were collected by the NORA historical ecology working group (extracted 27/01/2023). References to countries, provinces, regions, cities or offshore areas are mapped in Stechele et al. (2023). Designation of offshore regions is done according to Met OfficeSea Areas and Coastal Weather Stations referred to in the Shipping Forecast.

2.2 Larvae dispersal model

The LARVAE&CO model is an Individual-Based Model (IBM) that simulates egg and larval dispersal in the Eastern English Channel and the North Sea. It results from the coupling between a 3D hydrodynamic model and a Lagrangian particle-tracking module (Savina et al., 2010). This model, originally developed for sole and described in Lacroix et al. (2013) has been adapted for European flat oyster.

The 3D hydrodynamic NOS (North Sea) model, based on the COHERENS model (Luyten, 2006), has been implemented in the area between 48.5° N, 57° N, 4° W, 9° E. The model domain contains a 157 x 205 horizontal grid with a resolution of 5’ in longitude and 2.5’ in latitude (approximately 5 x 5 km) and 20 σ-coordinate vertical layers. Two open boundaries are located at the northern and western limit (at 4°W and 57°N) and the model includes daily river discharges of 14 rivers. The model is forced by weekly sea surface temperature (SST) data on a 20 × 20 km grid interpolated in space and time according to the model resolution (Bundesamt für Seeschifffahrt und Hydrographie, BSH, Germany) (Loewe, 2003) and by six-hourly surface wind and atmospheric pressure fields (provided by the Royal Meteorological Institute of Belgium based on the forecast data of the UK Met Office Global Atmospheric Model) (Hi Res, Walters et al., 2017). Details about the model implementation can be found in Savina et al. (2010) and Lacroix et al. (2013).

Larval trajectories were calculated online using the particle tracking model. The vertical diffusion was modelled by the random walk technique (Visser, 1997). Because in the North Sea vertical turbulent diffusion is considered to be the dominant horizontal diffusion mechanism (Christensen et al., 2007), explicit representation of horizontal diffusion was neglected. Specific details on the implementation can be found in Lacroix et al. (2013). Because current oyster bed distribution is unknown, connectivity between sites was not considered and only self-recruitment was modelled.

2.2.1 The individual-based model

In this implementation, only one pelagic stage was considered representing the pelagic larvae, from the swarming up to the settlement. The swarming period was identified to be temperature dependent just like the Pelagic Larval Duration (PLD). Larvae were passive and neutrally buoyant: no swimming behaviour or vertical migrations were considered. No natural mortality other than that due to loss by dispersal was included.

2.2.2 Parameterization

2.2.2.1 Spawning grounds

All the grid cells of the North Sea were potentially considered as spawning areas to allow to determine potential dispersal. To determine spawning grounds, a clustering method was applied based on the day where the temperature reached 15°C (corresponding to the first day of fertilisation), which was computed each year (over the period 2000-2009) for all grid cells. In the case where 15°C was not reached in a particular cell for a specific year, the day of the highest annual temperature was chosen. Numbers of the first spawning days were grouped by clustering method. Composition and the numbers of time clusters were determined with KMEANS algorithm with the package mclust in R (Scrucca et al., 2016) (Figure 1).

FIGURE 1

Figure 1 Release areas for European flat oyster. Each color represents an area with a similar spawning period (cluster). Spawning periods are computed each year (over the period 2000-2009) on each grid cell, and then aggregated using a clustering method. Dots indicate locations of historical oyster beds (Thurstan, 2023).

2.2.2.2 Release period

Eggs are fertilized within the females when temperature reaches 15°C (reproduction) within the cluster. Larvae are then released after a brooding period of ten days (De Mesel et al., 2018 and refs in). Release of larvae occurs over the course of a period of 10 days (De Mesel et al., 2018and refs in). Release intensity increases during the two first days and decreases during the eight following days. The timing of release (Julian day of the year) for the researched years, the release clusters and the release intensity is added in the Supplementary Materials Data Sheet 1.

2.2.2.3 Pelagic larval duration

Larval stage duration is a function of temperature according to the equation $L P D = α T^{β}$ where T is the temperature (in degrees Celsius). Published data on the duration of the larval stages, obtained from laboratory experiments and in situ measurements (De Mesel et al., 2018 and ref in) were used to estimate α and β parameters: $P L D (d a y s) = 1025.315 T^{- 1.56}$ .

2.2.3 Set of simulations

The model has been run over a period of 10 years (2000-2009) corresponding with the time scale of the North-Atlantic Oscillation; this to be representative for the broad range of hydrodynamic and environmental conditions in the area under consideration. The number of particles released in every grid cell is 632, which corresponds to a total of 9.153256 ×10⁶ particles released in the model domain. This number has been chosen based on a prior sensitivity study to arbitrarily ensure a good representation of the dispersal pattern.

2.3 DEB-IBM population model

2.3.1 General

The Dynamic Energy Budget (DEB) theory (Kooijman, 2010) provides a mechanistical framework that enables understanding the interaction between individuals, populations and ecosystems and the environment. In DEB-IBM population models (Martin et al., 2012) each individual of the population is represented as a unique entity by a set of state variables, unique DEB parameters (to account for inter-individual variability) and unique attributes. State variables are updated based on DEB theory to represent the individual change (growth, development, reproduction, and fitness) over the lifecycle of the individual. Population dynamics emerge from individual behaviour and interaction between individuals such as reproduction.

The state variables of the individual include energy reserve, $E$ (J), structural volume, $V$ (cm³), reproductive energy, $E_{R}$ (J), maturity level, $E_{H}$ (J), aging acceleration, $\ddot{q}$ (d^-2), age-related hazard rate, $\dot{h}$ (d^-1) and age, $a$ (d). The DEB model which is used to update the state variables is described by Stechele et al. (2022) and is an asj-model that accounts for metabolic acceleration from settlement onwards. Besides these DEB state variables, individuals are characterized by an individual-specific zoom factor which accounts for inter-individual variability in metabolic performance (Stavrakidis-Zachou et al., 2019). The structural length at settlement, $L_{s}$ (cm) and the structural length at the end of metamorphosis, $L_{j}$ (cm) are also individual-specific attributes that define the metamorphosis (metabolic acceleration in DEB terminology). The list of the individuals and their state variables, parameters and attributes represent the population.

2.3.2 DEB-IBM

Simulations start with 10 healthy adults. The state variables of the individuals are updated daily based on the environmental forcings. Mortality is checked daily, and dead individuals are removed from the populations. Mortality occurs when the individual’s energy reserve is empty $(E < = 0)$ , when individuals are preyed upon ( $R < \dot{p}$ ) or when accumulated tissue damage is too high for survival ( $R < \dot{h}$ ), with $R$ a random number between 0 and 1 and $\dot{p}$ the non-specific predation probability. Reproduction, larval release and swarming is triggered by water temperature (see 2.2.3). During reproduction, female eggs are fertilized and develop into larvae during the brooding phase. The release of larvae is distributed according to the release intensity which is given in Supplementary Materials Data Sheet 1. The number of larvae that are released $(n_{r e l e a s e d}),$ and are self-recruited $(n_{r e c r u i t e m e n t})$ by the population, are calculated accordingly:

n_{r e l e a s e d} = E_{R_{p o p}} n_{f e m a l e} E_{0}^{- 1} \frac{# L a r v a e_{r e l e a s e d}}{632}

n_{r e c r u i t e d} = n_{r e l e a s e d} S_{b g} S_{d r i f t} (R T_{m i n} + R T)

Where, $E_{R_{p o p}}$ is the sum of the $E_{R}$ for all mature individuals of the population, $n_{f e m a l e}$ the fraction of individuals that are females, $E_{0},$ the energy content of one egg (J), and $# L a r v a e_{r e l e a s e d}$ according to the distribution given in Supplementary Materials Data Sheet 1. Recruitment depends on the number of released larvae, the survival of the larvae through their pelagic drift phase ( $S_{d r i f t}$ ), and during the first year of life ( $S_{b g}$ , Walne, 1961). Survival during the first year of life (background survival) is integrated in the larval mortality rates to save on computation time. $R T_{m i n}$ represents the fraction of larvae that stay in the pixel of release, while $R T$ represents the fraction of larvae that drift away from location of release but return to the pixel of release due to current patterns. Larvae that drift away of the pixel of origin and arrive in another pixel (connectivity) are not considered. Temperature and food availability during the 20 days prior to settlement is averaged and is used to quantify the state variable $E$ of the recruited larvae. The age and length at settlement, $a_{s}$ (d) and $L_{s}$ (cm) were computed using the dget_ts function of the DEBtool. Other initial state variables of the recruited individuals include $[V = L_{s}^{3}; E_{R} = 0; E_{H} = E_{H_{s}}; q = 6.7 10^{- 10}; \dot{h} = 6.7 10^{10}$ ]. An individual-specific zoom factor is given randomly to each newly settled individual according to a normal distribution of $z \pm z_{σ}$ . The settling larvae that are integrated in the population behave as the other individuals and follow DEB rules during further development. When young individuals reach the end of the metamorphosis, $L_{j}$ is recorded individually as an attribute of the individual. $L_{s}$ and $L_{j}$ define the acceleration factor and individuals that are underfed during the early development will therefore grow slower and reach a smaller ultimate length.

The DEB-IBM model is linked to the larvae dispersal model (LDM) through the following parameters: (1) the timepoints at which reproduction occurs, (2) the timepoints at which release occurs, (3), the timepoints at which larvae arrive back in the pixel of release and (4) the fraction of released larvae that return to the pixel of origin (self-recruitment). For each location, population size, average fitness of all individuals (population fitness), larval release, and recruitment were recorded on a daily basis. The DEB-IBM parameters are presented in Supplementary Materials Data Sheet 2. The model assumptions are given in Table 1.

TABLE 1

Table 1 Model assumptions of the coupled DEB-IBM-LDM.

2.4 Suitability indicators, validation, initial conditions and sensitivity analysis

When envisioning restoration, one is generally interested in allocating restoration efforts to locations where population growth is fast and suitable substrate is present. Pixels with high (i) population increase, expressed as final population size after 10 years of growth are therefore considered suitable for restoration. Nevertheless, historical beds are not necessarily located at locations with high population increase. Other indicators for suitable locations for restoration include (ii) population fitness expressed as minimum average population fitness of all individuals of the population over a 10-year period (Stechele et al., 2023), (iii) population reproductivity expressed as the total larval release by the population over a 10-year period, and (iv) population self-recruitment expressed as the total number of larvae (over a 10-year period) that, after drift, arrive in the same location as where they were released. These suitability indicators are outcomes of the model and can be linked to the success of certain flat oyster populations. The values of these population suitability indicators are categorized as high, average, and low (quantification is given in the Results section). Pixels with a high or average suitability indicator value receive a suitability indicator score of 1, while pixels with low suitability indicator values receive a suitability indicator score of 0. Suitability indicator scores are summed to reach a total suitability score for each pixel. For larvae to settle, and populations to develop, (v) suitable substrate is an additional key requirement. Oyster populations generally do well on coarse substrates (including shells), stable sands, and stiff muds that provide a solid foundation. Therefore, soft muds and shifting sandy substrates are unsuitable (Héral and Deslous-Paoli, 1991; Smaal et al., 2017; Rodriguez-Perez et al., 2019; Hughes et al., 2023).

The probability of encountering a historical bed in a pixel with a certain suitability indicator value was used as a way to validate the results of this study. The validation was performed for each suitability indicator and for the total suitability score.

A sensitivity analysis of the model was conducted, involving the variation of several parameter values, including larval survival during drift $S_{d r i f t}$ , non-specific predation probability $\dot{p}$ , background survival $S_{b g}$ , and standard larval retention $R T_{m i n}$ . During the sensitivity analysis, population parameters were subjected to a positive and negative 20% variation. The primary objective was to observe the impact of these variations on the suitability indicator values, which serves as the main model outputs. Also, the impact of the initial conditions on the initial suitability indicator values was evaluated. The effect of the introduced adults was tested across a range of 10 to 10,000 individuals. Notably, only first-order impacts were examined, with each parameter being altered independently.

The sensitivity analysis and the impact of initial population size was carried out at two distinct locations since they are location dependent. Location 1 (55.1250° N, 8.4167° E) was identified as a site with low population increase, and after a 10-years development period, the population size was estimated to reach 306 individuals. In contrast, Location 2 (50.0833° N, 1.4167° E) supported higher population increases, and following a similar 10-year growth period, the population is expected to comprise 4647 individuals.

Matlab R2020a and its Mapping toolbox were used to perform simulations and generate plots.

3 Results

3.1 Suitability indicators

Population densities are expected to change over time, with population increases (Figure 2A) ranging from 0 to 5.5 × 10¹⁰ added individuals, which corresponds to an increased density of 0 to 1.1 × 10⁸ oysters/km² over a period of 10 years, depending on the location. Large population increases (>10⁶ individuals) are expected to occur in sheltered nearshore areas in the Wash (UK), around Isle of Wight (UK), on the north coast of Finistère (FR), Côte Fleurie (FR) all along the coast of Haute Normandie from Fécamp to Le Tréport (FR), around Dunkerque (FR) and at some locations in the Eastern and Western Frisian Islands (DE). An average increase (10⁴ – 10⁶ individuals) in population size is expected in the vicinity of areas with high population increase, in the Firth of Forth (UK), the coastal and nearshore Thames estuary (UK), Baie des Veys (FR), north coast of Nord-Pas-de-Calais-Picardie (FR), in the Scheldt estuary (BE, NL), around Noord Holland, and around German West and East Frisian islands (DE). Offshore locations with an average population increase are expected to be patchy in the English Channel but to be absent in the North Sea except for the area around Helgoland (DE). Expected population increases in the offshore North Sea area in general are low (< 10³ individuals).

FIGURE 2

Figure 2 Expected distribution of suitability indicators (based on flat oyster populations during a 10-year period 2000-2010). (A) Population increases (log scale). (B) population fitness (0-1). (C) population reproductivity (log scale). (D) population self-recruitment (log scale). Dots indicate locations of historical oyster beds (Thurstan, 2023).

The population fitness (Figure 2B), expressed as the minimum fitness (0 - 1, averaged over all individuals of the population) over a 10-year period, is expected to range between 0.10 and 0.63. The population fitness is expected to be high (> 0.5) in the nearshore areas of the English Channel including around South Devon (UK), Dorset (UK), Isle of Wight (UK), the north coast of Brittany (FR), the wider Bay of Saint Malo (FR), around the Channel Islands (UK) and the coastline of Haute Normandie and Picardie (FR). Population fitness is average (0.3 - 0.5) in most areas of the English Channel (Plymouth, Portland, Wight and Dover) and in some areas of the North Sea (Tyne, Humber and Fisher). Low population fitness (< 0.3) is expected to occur in large areas of the North Sea including the offshore areas of Thames, the German Bight, Forties, Dogger and Forth.

The populations reproductivity, expressed as the number of larvae released by the population over a 10-year period (Figure 2C), is expected to range between 1.1 × 10¹⁴ and 2.0 × 10²⁰. Total population reproductivity is uniform in the English Channel (Plymouth, Portland, Wight, Dover), the Southern Bight of the North Sea (Thames, Humber), and the German Bight. Larval production is high (> 10¹⁷) along all the coastlines, except for the Elbe estuary and in nearshore areas, while being average (10¹⁶ – 10¹⁷) in the offshore areas, except for the offshore areas of Forth, Forties, Dogger and Fisher (< 10¹⁶) The self-recruitment (Figure 2D), expressed as the number of recruited larvae over the 10-year period, ranges from 2.0 10² to 4.7 10¹². High self-recruitment (>10⁷) occurs in the Wash (UK), in the Thames estuary (UK), around Isle of Wight (UK), the Bay of Saint-Malo (FR), the Côte Fleurie (FR), all along the coast of Haute Normandie (FR), the north coast of Nord-Pas-de-Calais-Picardie (FR), all along the Belgian and Dutch coast and around the North and West Frisian Islands (DE). Average self-recruitment (10⁵-10⁷) is patchy in the western offshore English Channel (Plymouth, Portland), in almost all the offshore areas of the eastern English Channel, all over the Southern Bight of the North Sea, except for the territorial waters of Belgium and The Netherlands (with the exception of the coastal areas as mentioned above). Low self-recruitment (< 10⁵) is expected to occur along the North Sea coast of West-England (Tyne), the central and northern parts of the North Sea (Dogger, Forties and Forth), and off the coast of northern Denmark (Fisher). All datasets can be found in Supplementary Materials Data Sheet 3.

3.2 Suitability index score and validation

The total suitability indicator score (number of suitability indicators that are average or high) evaluates the suitability of locations for flat oyster populations. The flat oyster niche can be based on population dynamics (Figure 3A), or a combination of population dynamics and suitable sediment (Figure 3B).

FIGURE 3

Figure 3 Spatial variability of the total suitability indicator score based on population dynamics (A) and both sediment type and population dynamics (B). Dots indicate locations of historical oyster beds (Thurstan, 2023).

Restoration hotspots (total suitability indicator score = 4) are located in the Wash (UK), all along the coast of East England and the outer Thames estuary (UK), in the offshore mid Channel south of the Island of Wight (Wight), around the Isle of Wight (UK), all along the coast of Dorset and Devon (UK), all along the coast from Saint-Brieuc (FR) to Le Mont-Saint-Michel (FR), Côte Fleurie (FR) all along the coast of Haute Normandie from Fécamp to Le Tréport (FR), around Dunkerque (FR), in the Northern parts of the Scheldt Estuary (NL), off the coast of Noord-Holland (NL), in and around both the Western and Eastern Frisian Islands (DE) and around Helgoland (DE).

Suitability indicators scores can be linked to the occurrence of historical beds. Historical beds generally (>5% coverage) collocate with locations that have a high total suitability indicator score (Figure 4). There is a 10% chance of encountering a historical bed at a location with a total suitability indicator score of 3, while there is a 33% probability of encountering a historical bed at a location with a total suitability indicator score of 4 (Figure 4E). In locations with lower total suitability indicator scores, one has a lower probability of finding historical habitats. This both validates the results of this study and might explain why some historical habitats prospered at certain locations.

FIGURE 4

Figure 4 Validation of suitability indicator values against presence of historical oyster beds (probability of finding an oyster bed in a pixel with a certain suitability indicator value). Suitability indicators include population increase (A), population fitness (B), population reproductivity (C), population self-recruitment (D), and total suitability indicator score (E). The probability of finding oyster populations also depends on sediment type (F).

There is a probability of 6 to 33% to encounter historical beds in pixels with an average or high population increase. The probability of encountering a historical bed increases when the expected population increases. Similarly, the probability of encountering historical beds increases for pixels with high population fitness. In pixels with a minimum population fitness > 0.6, one has a 25% probability of encountering historical beds. Historical beds generally occurred in location where larval production is estimated to be average or high (>10¹⁶) and where population self-recruitment is average or high (10⁵). In relation to substrate type, historical habitats coincide most with areas composed of rocks, boulders, coarse substrate, mixed sediment, mud and sand, in order of priority, and most likely aligns with the stability of the substrate as well.

3.3 Sensitivity analysis and the effect of initial population size

The model results are influenced by forcings, model parameters and initial conditions. The influence of population parameter values on the suitability indicators is given in the sensitivity analysis (Table 2). Population fitness is generally unsensitive to population parameter variations. The other suitability indicators do are sensitive to variations in population parameter values.

TABLE 2

Table 2 Sensitivity analysis for the DEB-IBM population parameters for a location with low population increase (Location 1; 55.125° N, 8.4167° E) and average population increase (Location 2; 50.0833° N, 1.4167° E).

Variations in $S_{d r i f t}$ result in a low to moderate impact on the population increase, reproduction, and self-recruitment. These suitability indicators are highly sensitive to variations in $\dot{p}$ and $S_{b g}$ and unsensitive to variations in $R T_{m i n}$ . Model outcomes generally become more sensitive to variations in parameter values when the environmental conditions are more suitable.

The population size at the start of the simulation affects the model outcomes. The positive effects of additional broodstock introduction becomes more important when environmental conditions are suitable. For a location where low population increase is expected (Location 1), increasing the initial population size from 10¹ individuals to 10⁴ resulted in a doubling of the final population size (from 306 to 595 final population size). In a more suitable location (Location 2), on the other hand, population size increased with a factor 42 (from 4647 individuals to 197104 individuals). The impact of initial population is therefore location dependent.

4 Discussion

4.1 General

The performance of organisms is strongly connected to their physical environment and their ability to cope with the fluctuations in their environment. Knowledge about the ecological niche of a species has been a key element in conservation and restoration efforts (Peterson, 2003; Peterson and Robins, 2003).

General trends in site suitability for flat oysters, not for flat oyster populations, are given in Stechele et al. (2023) and include; (1) coastal areas are rich in nutrients and support fast growth and high fitness, (2) low salinities due to river runoff are unsuitable, (3) offshore locations have more constant temperature profiles but lower food availability, (4) environments under Atlantic influence are characterized by low nutrients and therefore do not support flat oyster populations, (5) locations at higher latitudes in the northern North Sea, have low temperatures and do not reach the spawning threshold.

By extending the species model to a population model, the results of Stechele et al. (2023) are finetuned. The new trends that emerge restrict suitability in comparison to the Stechele et al. (2023) results. Spatial variability in site suitability for restoration depends on the environmental conditions at the site, including the food availability, salinity, and temperature dynamics as well as substrate type and hydrodynamics. Analyzing these factors on a site-by-site basis is needed to understand the reason why the location is suitable for restoration. However, the site suitability indicators (population increase, population fitness, reproductivity and self-recruitment) provide more understanding.

The high suitability for flat oysters (Stechele et al., 2023) all along the French coast of Brittany and Normandy translates into a patchy pattern of suitability for restoration. The coastline between Saint-Malo and the Mont Saint Michel for example, is estimated to be more suitable than the Côte des Isles in the west of Normandy. The reason for this can be found in the suitability indicators, and in this case relates to lower reproductivity and lower self-recruitment.

The use of ecological niche models and species distribution models allows one to explore the patterns behind observed species distribution. A comprehensive understanding of species distributions and site suitability is essential for successful habitat restoration initiatives which have grown in prevalence recently (Sillero et al., 2021). Correlative approaches are typically employed to evaluate species distribution; however, correlative niche models lack the capacity to explain the processes behind the observations and are thus limited in various applications, such as site selection for habitat creation (habitat creation involves the creation of ecosystem at locations where these systems previously did not occur), evaluating the effect of climate change on niche distribution, and understanding migration patterns of invasive species. In general, these models lack physiological knowledge as to why certain factors are linked to the presence or absence of a species (Kearney and Porter, 2009). Mechanistic niche models, as developed and applied in this study, are of particular importance when one envisions understanding flat oyster habitats distribution because (1) precise records of past bed locations are often incorrect, (2) the historical environment has changed significantly and will continue to change and (3) there is an increasing interest in flat oyster nature-inclusive-designs with the aim of restoring flat oyster habitat functions on locations where protection measures have been in place (e.g., offshore wind farms, MPAs, offshore infrastructures).

4.1.1 Advantages of the methodology

This work presents a novel method of evaluating site suitability using a realistic niche model that is supported by a generic and mechanistic approach. The generic nature of DEB implies that the methodology proposed here is applicable to other species (Kooijman, 2010; Kearney, 2012). Species specific parameters can be found in the Add-My-Pet collection, and by applying these, one can use the same modelling approach for other sedentary species. The mechanistic nature of the DEB model, which forms the basis of this current ecological niche model, has the capacity to explain the behaviour of individuals and population in variable environmental conditions. Simulating state variables of individuals in a population, describes processes such as growth, starvation, mortality and successful or failed reproduction events. By including inter-individual variability and individual interactions one surpasses the individual level. Mechanistic population models become relevant in explaining population traits such as population dynamics, population density, reproductivity or self-recruitment as well as quantifying regulating ecosystem services (water quality regulation and nutrient removal) and evaluating supporting services (habitat provision, biodiversity, sediment biogeochemistry alterations and food web structure) (Rullens et al., 2019).

The proposed methodology additionally has the benefit of being straightforward to comprehend, modify, or enhance. DEB extensions to include the effects of hypoxia (Thomas et al., 2019), filtration inhibition by particulate matter (Kooijman, 2010), or gonadal development (Bourlès et al., 2009) would increase the explanatory power of the model, but these DEB extensions have not yet been parameterized for flat oysters. The DEB model for flat oyster individuals has been parameterized and validated (Stechele et al., 2022).

As done in this study, these mechanistical population models can easily be validated using presence/absence datasets, in contrast to correlative niche models which often stay unvalidated (Theuerkauf and Lipcius, 2016).

4.1.2 Limitations

Although the methodology is applicable on all time and space scales, the application discussed in this work is limited by both scales. The applied timeline (10-years) accounts for the North-Atlantic Oscillation’s typical variability. It does not account for rare environmental events that might have an implication on flat oyster populations (e.g., the severe winter of 1962/1963 caused a reduction of the flat oyster stock by 97% in the Eastern Scheldt estuary, rare storm surges or current shifts have implications on water temperature or plankton distribution) (Drinkwaard, 1998; Smaal et al., 2022). Spatially, the model is applied at mega (basin) and macro scale (km). Suitability at meso (m) and micro (>m) scale are important to consider and are known to impact model recruitment parameters ( $R T$ and $R T_{m i n}$ ) (Hofstede et al., 2023). Applying the model to finer resolution data is essential when evaluating site suitability characterized by steep gradients such as estuaries or dynamic coastal environments (Bergquist et al., 2006) or even at smaller scales such as artificial habitats (Bartol and Mann, 1997).

Another important limitation of this study is the valuation of important parameters. For example, mortality of individuals is assumed when their energy reserve (fitness) is empty ( $E = 0$ ). Oysters are known to re-allocate reproductive energy to cover maintenance needs during starvation periods (not included in this model) (Flye-Sainte-Marie et al., 2009; Liu et al., 2010). On the other hand, weak individuals might be more susceptible to disturbances (such as temporary burial, disease, exposition to toxic algae or pollution) and mortality might occur at fitness levels higher than zero ( $E > 0$ ) (Mcdonald et al., 2023). Valuation of the recruitment parameter, $R T_{m i n}$ depends on meso- and microscale relief of the site and the substrate and is therefore difficult to assess, although it can influence the suitability indicator population self-recruitment. Another important parameter which is difficult to value is the fraction of individuals which are female ( $n_{f e m a l e}$ ), a parameter that impacts larval release (impacts suitability indicator population reproductivity). The flat oyster is an asynchronous hermaphrodite organism which changes sexuality over the year. The mechanism of sex change is poorly understood (da Silva et al., 2009). To overcome the limitations related to the valuation of important parameters, suitability indicators have been developed in this study. These suitability indicators are closely linked to parameters that are difficult to quantify. Locations that support good population fitness (fitness > 0.4), are more likely to harbor successful oyster beds restoration, since these beds are well fed and therefore more resilient to environmental disturbances. Similar reasonings can be followed in case of the other suitability indicators. The sensitivity analysis clearly shows that the valuation of the population parameters is important when using quantitative model results. The model outcome results are highly sensitive to variations in the parameters $S_{d r i f t}$ , $\dot{p}$ and $S_{b g}$ , indicating that more research is needed to quantify and understand survival and predation characteristics of sites.

The relation between the occurrence of historical habitats and Folk 5 substrate types (mud, sand, mixed sediment, coarse substrate, rocks, and boulders) generally aligns with the degree of stability that is expected from these beds. Whilst this is clear for the hard substrate types, more variables are into play for the soft substrates. In Folk 5, the sand class is most narrowly defined as a sediment with less than 10% mud and no gravel fraction, hence it is the most mobile, and indeed least suitable. The mud class can have up to 90% mud and can vary largely in composition (muddy sand to clay). Importantly, one should consider the substrate beneath the active layer, that is often very thin (mm to cm) (van der Spek et al., 2022), and overtopping more compact, and often older geological layers that would provide the necessary stability and firmness. Depending on the thickness of the active layer, time dependent seabed samples represent approximatively these more compact layers. Hitherto, this is not mappable at larger scales. Mechanistically accounting for seabed stability is difficult since it implies accounting for several hydrodynamic factors such as tidal and residual currents, depth, and storm surges, together orchestrating the dynamicity of the seafloor geomorphology and the substrates composed of variable grain sizes. Seabed mobility can be approximated using proxies such as depth, bed shear stress or mud content in sandy sediment (Damen et al., 2018). The percentage of mud in a sandy sediment may be indicative of relatively stable conditions, if consolidation can occur, e.g., where fine silt and clay particles can settle (Borsje et al., 2013; Damen et al., 2018). Flat oysters appear to have a preference for consolidated sediment with a relatively high mud content. However, seabeds with extremely high (>50%) unconsolidated mud is unsuitable (Houziaux et al., 2008; Houziaux et al., 2011).

However, the regionally available EMODnet products result from standardization and harmonization of many datasets, using a range of methodologies and often based on sampling data which cover a wide time span (Kaskela et al., 2019). Significant improvements are expected with upcoming EMODnet products which are based on more coherent datasets and more advanced methodologies. Progress is also made in a better parameterization of the geological substrate, and sediment thickness (e.g., Hademenos et al., 2019). Incorporating better sediment parameterization in seabed mobility studies will further result in more realistic suitability predictions. Considering correlative mismatches between areas predicted as suitable and the presence of historical habitats (false positives), one should consider the reliability of the substrate/seafloor terrain variables used in this study, other missing variables, and the correctness of historical oyster bed distribution data, which were mainly drafted when overfishing was already ongoing for more than a century and the disappearance of beds was had being recorded.

Therefore, it must be noted that the EMODnet substrate dataset used in this study may not be sufficiently representative of such complex and important dynamic factors. It must also be noted that there exists a paucity of datasets with large-scale geographical coverage as applied in this study.

Observing the early life stages of oysters directly can be difficult at the North Sea scale, as they are tiny. Therefore, biophysical models are used to assess their dispersal, which implies the formulation of hypotheses about spawning periods and pelagic duration. However, the validity of these hypotheses can be discussed, especially in areas where populations are not naturally present, such as the northern part of the domain where, in the model, spawning occurs below the temperature threshold of 15°C, Nevertheless, these areas also have low suitability indicators making the hypothesis less critical (Figure 3).

The size of the grid choice also plays an important role in measuring the self-recruitment of oyster populations. Indeed, the smaller the grid size, the fewer particles will drift and remain in place. This can lead to an underestimation of self-recruitment if the grid size is too small, as particles can drift outside the grid and contribute to recruitment in another area. However, a too large grid size could overestimate self-recruitment. The grid size used (5 km x 5 km) is a compromise that balances the need to capture fine-scale details with the uncertainty in particle movement. While such a grid may not be representative of the specific restoration sites, it still provides valuable insights into self-recruitment dynamics on a regional scale.

This study, which focuses on self-recruitment only, does not consider the potential connectivity with outside areas, be it natural, cultivated or restored populations. For instance, an area might be highly favourable for oyster restoration, even with a low level of self-recruitment provided that it receives spat from outside areas, and potentially it could even be a source of spat for other areas. Taking into account all potential connections would require a thorough knowledge of the existing population (location, density, population structure) and a network strategy for designing restoration sites across the North Sea. Connectivity was not considered in this work because nearshore and offshore population are heavily degraded. Locations and densities of populations in the nearshore and offshore environment are unknown. In addition, the aim of restoring beds is to generate fully functional, self-sustaining beds. Dependence on nearby populations would make restoration sites less resilient.

An important model assumption (Table 1) is that population size does not influence the environment of the site. Oyster populations generally influence the environment by regulating food levels (Akinyemi et al., 2021; Albentosa et al., 2023), providing suitable substrate for settlement (Rodriguez-Perez et al., 2019), increasing complexity of the seabed which improves larval retention (Andutta et al., 2012) or attracting biodiversity including oyster predators. These relationships with the environment define carrying capacity or maximum population size. Although maximum population size is an important characteristic in restoration and site selection studies, it is not represented in the results of this study, and the results of this study merely relate to the rate at which populations can be restored. Large population increases (>10⁶) as estimated by this study, correspond to small population densities (> 0.04 oysters m^-²), assuming even distribution of oysters over the pixels, after 10 years. These densities do not represent oyster beds or reefs, but merely represent sporadic individual oysters which only minimally impact on food levels, substrate availability or biodiversity (zu Ermgassen et al., 2021).

4.2 Recommendations to restoration practitioners and outlook

As indicated in this study, the model provides several outputs that can serve as indicators for a scientifically informed estimation of oyster habitat restoration success. A rapid population increase is a desirable factor for restoration initiatives and is a keystone indicator of success (Glover et al., 2021). Here, high population increases are estimated to occur within sheltered and nearshore areas where reproductivity and larval retention are high. Offshore sites (especially in the North Sea) are generally characterized by low larval retention. We do not expect to see high population increases at these sites. Supplementing restoration efforts with high numbers of broodstock individuals seems self-evident at offshore sites, but hydrodynamics disperse larvae rapidly, with limited expected influence on success. Management strategies that would have more influence in the offshore environment is reducing predation, increasing complexity, to keep larvae on site, or applying large numbers of spat on shell. Scaling up restoration efforts will boost population productivity and self-recruitment.

To increase success of flat oyster habitat restoration in the offshore environment, we suggest the implementation of a basin-wide coordinated restoration effort that promotes the connectivity between natural oyster beds, restoration sites, oyster nature-inclusive-design adaptation to offshore infrastructure and aquaculture sites. Allocating these supportive activities, to locations where high self-recruitment (or high connectivity) is expected, would benefit local beds.

The method presented in this work provides the basis to quantitatively understand why certain locations are more suitable for flat oyster habitat restoration than others and can be used to scientifically gauge environmental restoration efforts. In additions, the method represents a powerful tool for investigating multiple important research gaps including the potential impacts of climate change on oyster restoration efforts, the effect of restoration efforts on the ecosystems functioning (Kotta et al., 2023) or evaluate management strategies on population dynamics.

5 Conclusion

With the growing trend of flat oyster habitat restoration, creation, and increased importance of nature inclusive designs that integrate flat oyster habitats with offshore marine infrastructure, correct site selection through suitability mapping becomes an important prerequisite of success. This work demonstrates how combining an individual-based population model (DEB-IBM) and a larval dispersal model can be used to create a niche model for flat oyster habitat restoration. Spatial variability in suitability indicators; population increase, minimum population fitness, population reproduction and population self-recruitment define the total suitability. Through the application of this model to the English Channel and the North Sea, hotspots for restoration were identified, and reasons for restoration success or failure are identified through valuation of the suitability indicators. Additionally, the model can be used to explain why certain sites are suitable for flat oyster habitats and why certain historical habitats existed in certain locations. Furthermore, the model can be used to measure the effects of management strategies and the ecosystem services provided by habitats.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics statement

The manuscript presents research on animals that do not require ethical approval for their study.

Author contributions

BS, LB, GL, LD and NN conceived the ideas and designed methodology. BS, LB, GL built the databases. BS, LB, GL, LD, GG analyzed the data. BS, LB, GL, VL and GG wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the Research Foundation–Flanders (FWO) through an SB PhD fellowship granted to BS [project number 1S84619N]. This project has received funding from the European Union’s Horizon 2020 Research and Innovation Programme under Grant Agreement no 862915. Views and opinions are those of authors only and do not necessarily reflect those of the European Union. Neither the European Union nor the Granting Authority can be held responsible for them.

Acknowledgments

This study has been conducted using E.U. Copernicus Marine Service Information; https://doi.org/10.48670/moi-00054 and https://doi.org/10.48670/moi-00058.

Conflict of interest

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

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

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

References

Akinyemi M. I., Da-Costa F., Filgueira R., Galimany E., Ibarrola I., Albentosa M. (2021). Scope for growth of Mar Menor flat oysters under eutrophic conditions: Implications for future bioremediation actions. [Conference presentation]. NORA 4.

Google Scholar

Albentosa M., Akinyemi M. I., Vera M., Ibarrola I., Filgueira R., Galimany E., et al. (2023). Recovery of eutrophized marine ecosystems using the European flat oyster, Ostrea edulis. Aquat. Conserv.: Mar. Freshw. Ecosys. 33 (7), 645–660. doi: 10.1002/aqc.3926