# Ecological niche modelling of five gorgonian species within the shallow rocky habitat of the french mediterranean coast

Sylvain Blouet, Lorenzo Bramanti, Katell Guizien

► To cite this version:

Sylvain Blouet, Lorenzo Bramanti, Katell Guizien. Ecological niche modelling of five gorgonian species within the shallow rocky habitat of the french mediterranean coast. Vie et Milieu / Life & Environment, 2024, 7 (1-2), pp.1-43. 10.57890/VIEMILIEU/2024.74-1/2:1-44. hal-04637260

HAL Id: hal-04637260

https://hal.science/hal-04637260v1

Submitted on 10 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ECOLOGICAL NICHE MODELLING OF FIVE
GORGONIAN SPECIES WITHIN THE SHALLOW
ROCKY HABITAT OF THE FRENCH
MEDITERRANEAN COAST

S. BLOUET<fnref n="1" />,<fnref n="2" />*, L. BRAMANTI<fnref n="1" />, K. GUIZIEN<fnref n="1" />

GORGONIANS
MAXENT
MARINE PROTECTED AREA
DISTRIBUTION
HABITAT SUITABILITY
MEDITERRANEAN

Abstract.- Engineer species, such as gorgonians, provide several ecosystem services and play a significant role in the maintenance of biodiversity. We tested an ecological niche modelling approach using a unique dataset of spatialized inventories on a regular grid (< 800 m) replicated in 2013 and 2020 along 450 km of coastline for five gorgonian species commonly found in the Mediterranean shallow habitats (10-50 m deep). Ten non-correlated environmental predictors derived from the most advanced geomorphological and hydrological data were used to assess the ecological niche of the five species using maximum entropy (Maxent). The model confirmed that different combinations of depth, rugosity, sea surface and bottom temperature, bottom flow speed and turbidity explain the different gorgonian's distribution. The reduction of the number of presence observations did not alter the quality of the prediction as long as the observation points were spread over the entire variability range of predictors. The latter can be achieved by only including presence observations from highly protected zones of the marine protected areas of the region. Our results provide a greater understanding of the factors shaping the distribution of five gorgonian species commonly found in Mediterranean shallow areas.

INTRODUCTION

The maintenance of the ecosystem services
provided by the ocean and that contribute to hu-
man welfare is one of the 17 sustainable develop-
ment goals identified by the General Assembly of
the United Nations in 2015 (https://sdgs.un.org/
fr/goals). Such an objective is encompassed in
the European Marine Spatial Planning framework
which aims to take into account the interactions
between different human activities at sea, in-
cluding biodiversity protection (DIRECTIVE
2014/89/UE). One-step in marine spatial planning
is mapping biodiversity and the ecosystem services
it delivers, together with understanding its drivers,
in order to be able to forecast its evolution under

different development scenarios (Bailey 2010).

During the last 20 years, in the Mediter-
ranean Sea, several research programs focused on
obtaining spatial information on marine habitats
and the distribution range of marine species in or-
der to begin implementing marine spatial planning
(Gerovasileiou et al. 2019). In this region, marine
species inventories, however, have generally been
restricted to a limited number of locations, often
in Marine Protected Areas (MPAs), although re-
cent citizen science projects attempted to expand
this spatial scale (Bramanti et al. 2011, Ponti et al.
2019). Ecological niche modelling (ENM) can be
a powerful complementary tool to refine species

<footnote n="1">¹ CNRS-Sorbonne Université, Laboratoire d'Ecogéochimie des Environnements Benthiques, LECOB, Observatoire Océanologique de Banyuls Sur Mer, 66650 Banyuls sur Mer, France</footnote>

<footnote n="2">² Ville d'Agde, Aire marine protégée de la côte agathoise, 34300 Agde, France</footnote>

<footnote>* Corresponding author: sylvain.blouet@ville-agde.fr</footnote>
habitat mapping. It currently relies on correlative approaches, which are based on the assumption that the suitable environment of a species controls its range and allows one to project species distribution range in the geographical space using environmental predictors (Sillero et al. 2021). Early ENMs used statistical modelling such as Generalised Linear Modelling (GLM, McCullagh et al. 1983) relating presence-absence to environmental descriptors. Machine learning models that use presence-only data (Maxent, Phillips et al. 2006, Random Forest, Breiman 2001), together with satellite imaging boosted ecological niche modelling applications at large spatial scale. The performance of various ENM algorithms have been systematically compared, showing species-related performance and suggesting no single algorithm outcompete others strengths and limitations (reviewed in Guisan et al. 2017). An important conclusion about the performance of ENM is their sensitivity to the number of environmental predictors used (Sainz-Villegas et al. 2022). The ensemble averaging of multiple ENM outputs has been proposed to improve species spatial distribution mapping (Achu et al. 2021). Still, this improvement blurs the explanation link of the environmental predictors. Among ENMs, the machine-learning model Maxent is one of the most tested. It displays ecologically meaningful responses to environmental predictors and is less data (presence-only) demanding than GLM, its statistical counterpart that requires more extended datasets with presence and absence (Elith et al. 2006, Wisz et al. 2008, Sainz-Villegas et al. 2022). Specifically, using presence-only data avoids the assumption that the absence of a species means habitat unsuitability, when in fact it could be the result of either disturbance or ineffective dispersal (Jarnevich et al. 2015). The most commonly reported sampling biases are geographic biases - such as partial distributions and environmental biases - leading to either unreliable results or to niche estimates with underrepresented subsets of the population environmental space (Sillero et al. 2021). The pressures on the habitats can also bias model construction. Indeed, the absence of a species can result from anthropogenic pressures irrespective of the characteristics of the habitat to support this species, that should be disregarded in ENM appli-

cations. Highly protected zones within MPAs might therefore be preferred locations for presence observations in order to limit confusion between predictors of the ecological niche and habitats degraded by direct anthropogenic pressures. In the present study, we aimed to develop an ecological niche model for the most common gorgonian species in the NW Mediterranean, using an extensive dataset of presence observations covering a wide range of environmental conditions in the NW Mediterranean Sea, inside and outside of MPAs.

Gorgonians (octocorals with an internal skeletal axis, McFadden et al. 2022) are among the most emblematic species of the Mediterranean Sea, playing an essential role in shaping and structuring the biological diversity on the rocky temperate marine habitats (Gili & Coma 1998, Otero et al. 2017). Among the Mediterranean gorgonian species, five are common in the rocky communities of shallow (0-50 m) coastal waters (Weinberg 1979a): *Eunicella singularis* (Esper, 1791), *Leptogorgia sarmentosa* (Esper, 1789), *Eunicella cavolini* (Koch, 1887), *Paramuricea clavata* (Risso, 1827) and *Corallium rubrum* (Linnaeus, 1758). The five species have a wide Mediterranean Sea distribution but also extend (except *E. cavolini*, López-González, 1993) to the Atlantic Ocean on the coasts of Morocco (Zibrowius et al. 1984, Grasshoff, 1977) and Portugal (Cúrdia et al. 2013, Boavida et al. 2016).

In the Mediterranean Sea, they have mainly been studied in the northwestern basin (Carpine & Grasshoff 1975, Cocito et al. 2002, Gori et al. 2011a, Di Camillo et al. 2018) although their distributions extend to the Tunisian and Sicilian coasts of the Ionian basin (Canessa et al. 2023), in the Aegean Sea (Salomidi et al. 2009, Vafidis, 2009, Çinar et al. 2018) and the Sea of Marmara (Topçu & Öztürk 2015). No populations were found in the easternmost part of the Mediterranean basin but this absence could be due to a lack of data (Otero et al. 2017). *Paramuricea clavata*, *C. rubrum* and *E. cavolini* have a wide bathymetric distribution from 10 m to more than 1000 m water depth (Grinyó et al. 2016, Knittweis et al. 2016, Carugati et al. 2022).

In this study, we analysed the spatial distribution of *E. singularis*, *L. sarmentosa*, *E. cavolini*, *P. clavata* and *C. rubrum* in their shallow bathymetry.
tric range (10-50 m) in the eight major rocky ha-
bitat spots of the French Mediterranean coasts
(spanning 3.5 degrees of longitude, from Toulon to
Cerbère), using spatialized inventories performed
on a regular grid (< 800 m). The ecological niche
and distribution area of suitable habitats were mo-
delled with Maxent software to identify the main
environmental predictors that influence the spa-
tial distribution of the five species. We tested the
sensitivity of the fitted model to a reduced num-
ber of occurrences and their geographic position,
and specifically assessed whether sampling limited
to the five highly protected areas of MPAs in the
study area (which can be considered exempt from
the main direct anthropogenic pressures impac-
ting gorgonians such as anchoring, nets, scuba
diving, etc.) would affect the model predictions.

MATERIALS AND METHODS

Study area

The study area extends along 450 km of the French Medi-
terranean coastline, from Toulon to Cerbère, and from the
coast down to the 50 m isobath. Along this coastline, the
seabed mainly consists of a soft bottom with a few small
patches of rocky habitat between 76 (VLR) and 3123 hec-
tares (AGM). A set of eight sites (National Park of Port Cros:
PNPC, Parc Marin de la Côte Bleue: PMCB, Aigues-Mortes:
AGM, Aresquiers: ARES, Agde: AGD, Valras: VLR, Leu-
cate: LEU, Côte Vermeille: CVM, Fig. 1) within this frag-
mented rocky habitat were investigated to assess the niche
preferences of the five species: *Eunicella singularis*, *Eu-
nicella cavolini*, *Paramuricea clavata*, *Leptogorgia sarmen-
tosa* and *Corallium rubrum*. Of the eight sites, five include a
high protection zone (according to the classification of Horta
e Costa *et al*. 2016: Fully protected) with different surface
areas (PNPC: 1232 ha, AGM:100 ha, AGD: 310 ha, CVM:
69 ha, PMCB: 289 ha, Fig. S1, Supplementary Material).

Data collection

*Species distribution:* The population density of the five species was assessed between 2013 and 2020 in 696 evenly spaced geo-referenced locations covering the hard bottom habitat in all sites except in PNPC where locations were selected based on the known presence of the species *P. clavata* (Guizien *et al.* 2022, https://cardobs.mnhn.fr/, Fig. 1 and Fig. S1, Supplementary Material). The spacing between sampling points varied from 100 m to 800 m, depending on the bathymetrical steepness of each zone. Using a regular grid in order to avoid the bias of preferential sampling enabled us to estimate that no spatial auto-correlation was observed (Legendre *et al.* 2002, Merckx *et al.* 2011). For the five species, a strong nugget effect (Carrasco 2010) at eight meters

distance was observed on empirical variograms (Fig. S2,
Supplementary Material). In each georeferenced location, all
individuals taller than 2 cm (i.e., older than one year) were
counted either by direct visual census of autonomous scuba
divers or on photographs in four to nine quadrats of 1 square
meter. These surface areas wer
e larger than the minimal area of 2 square meters establi-
shed for shallow water Octocorallia communities (Weinberg
1978). The georeferenced database of the five species pre-
sence and absence was implemented in a GIS software (©
QGIS).

Habitat suitability model

The Maxent model: Maxent is a machine learning software recognized as the most efficient method in species distribution modelling when presence-only species records are available (Elith et al. 2006, Hernandez et al. 2006, Hernandez et al. 2008, Bargain et al. 2018). Natively, Maxent estimates Receiving Operator Curve, more explicitly called Relative Occurrence Rate (ROR, Merow et al. 2013) which is defined as the probability for an observed individual to be found in a set of environmental conditions by maximising its entropy probability distribution (i.e., the most spread out or closest to the uniform, Phillips et al. 2006). Hence, Maxent is a niche or habitat suitability model whose predictions, although represented on a spatial grid, relate to environmental predictors only, disregarding spatial location arrangements except if geographical predictors are explicitly included in the set of environmental descriptors (Elith et al. 2011). The Maxent version 3.4.4 was used. Setting up the Maxent model requires defining a background grid on which environmental predictors are known, features, regularization multipliers, sampling biases, output types and diagnostics outputs (see Supplementary Material for details). The quality of Maxent fits is primarily assessed by the Area Under the ROR Curve (AUC). The AUC ranges from 0 (fit worse than random) to 1 (perfect fit). Below 0.5 the prediction is not better than at random, and above 0.7 the prediction is reliable (Zweig & Campbell 1993).

Data on marine habitats were retrieved from CARTHAM
website (https://geo.data.gouv.fr) and the background grid
covered the only hard bottom habitat where the five octocoral
species exclusively dwell. Gridding the hard bottom habitat
in the eight sites resulted in a set of 8942 georeferenced 100
m x 100 m cells (background grid), according to availability
of environmental predictors data (Fig. S1). Overlaying the
8942 background grid cells with the 696 sampling sites (men-
tioned above) resulted in 444 grid cells where both environ-
mental predictors were defined and the five species presence/
absence was assessed. Across the 444 grid cells, *Corallium*
*rubrum* was present in 16 grid cells, *P. clavata* in 18 grid
cells, *E. cavolini* in 31 grid cells, *L. sarmentosa* in 62 grid
cells and *E. singularis* in 255 grid cells. The implementation
of habitats and masks was carried out in © QGIS software.
The environmental predictors were selected within all avai-
lable geomorphological and hydrological potential predic-
tors, to which were added latitude and longitude as explicit
<figure><img src="image_2.png" /><figcaption>Fig. 1.- Map displaying the number of sampling locations in each of the eight sites of the NW Mediterranean Sea. Black rectangles highlight the sampled sites. The diameter of the orange circles is proportional to the number of sampling points at each site.</figcaption></figure>

predictors in some modelling test.

Geomorphological predictors were calculated from raw
bathymetrical dataset that combined three-dimensional point
seedlings from Airborne Laser Imaging Detection and Rang-
ging System (LIDAR) for water depth less than 15 m (LIT-
TO3D, www diffusión shom fr) and multibeam echosounder
for water depth larger than 15 m, except in PNPC where LI-
DAR data extended down to 40 m water depth. LIDAR data
was collected between 2014-2015, except in the CVM zone
where it was obtained in 2011. Multibeam echosounder data
was assembled from multiple collection periods between
2005 and 2018 and from different operators (Seaview for
AGM, ARES and VLR, CEFREM for CVM, Semantic TS
for AGD, Copethec and Mesuris for PMCB). Density of
LIDAR data decreased from 10 to 0.1 points per square
meter, when water depth increased. Density of multibeam
echosounder data decreased from 25 to 4 points per square
meter. Seven geomorphological metrics commonly used for
habitat predictive modelling were computed for each cell of
the background grid. These seven metrics describing the sea-
floor topography were: bottom depth (unit: m), slope (unit:
angle degrees), roughness (unit: m), rugosity (no unit), ter-
rain ruggedness index (TRI, unit: m), eastness (no unit) and
northness (no unit) orientation (Wilson et al. 2007, see Sup-
plementary Material for details). These metrics were calcu-
lated using all raw bathymetrical data within a 20 m x 20 m
slab located in the centre of each cell of the background grid.
Hydrological predictors included turbidity and tempe-

rature, which were measured at the sea surface by satel-
lite imaging (Aquamodis - level-3 global browser: https://
oceancolor.gsfc.nasa.gov/l3/), and bottom temperature and
flow speed, simulated by a coastal ocean circulation mo-
del (Briton et al. 2018). Satellite observations came from
the Moderate Resolution Imaging Spectroradiometer (MO-
DIS) on board the AQUA satellite. Data was acquired in 36
spectral bands at three resolutions (250, 500 and 1000 m)
and covered all the globe in one to two days. Sea Surface
Temperature (SST, unit: °C, 4 µm night time) and turbidity
(unit: m-1, Diffuse attenuation coefficient at 490 nm) month-
ly averages at 4 km spatial resolution were extracted between
2003 and 2018 for Turbidity and 2006 and 2018 for SST
over the extent of the background grid. For each predictor
(SST, Turbidity), the pluri-annual mean was calculated and
interpolated at each background grid point using the gridda-
ta Matlab function with the default linear option. Although
Chlorophyll-a concentration is one of the MODIS products
(Chla, unit: mg.m-3, OCx Algorithm), it was not conside-
red because data was missing in the PMCB site. However,
Chlorophyll-a concentration significantly correlated with tur-
bidity in all other sites (Spearman determination coefficient
R² = 0.68, p < 0.0001). Bottom temperature and flow speed
came from hourly outputs of a high-resolution SYMPHONIE
2015 regional circulation simulation, which ran from Janua-
ry 2010 to December 2012 (80 to 300 m horizontal resolu-
tion in the eight sites, Briton et al. 2018, Vissenaekens et al.
2023). Monthly statistics of flow speed (tenth and ninetieth
percentiles, unit: m.s-1) and temperature (tenth and ninetieth percentiles, unit: °C) hourly values at the bottom were computed. Those statistics were then interpolated in each cell of the background grid. As the horizontal resolution at the coast of the ocean circulation model ranged from 80 m to 350 m and water depth was thresholded to 2 m (Briton et al. 2018), the area where ocean circulation was computed was not perfectly overlapping the background grid depicting the rocky habitat near the coast. In this shadow area, ocean variables are undefined and interpolating values on the background grid thus required to define temperature and flow speed values at the coastline. Flow speed values were set to 0 at the coastline, meaning total dissipation of flow energy at the coast. Sea bottom temperature values were set to their average along the 2 m isobath of the ocean circulation model. Inter-annual bottom flow speed and temperature minimum (Umin and Tmin) and maximum (Umax and Tmax) at each cell of the background grid were then calculated, taking the minimum and maximum values of monthly tenth and ninetieth percentiles, respectively, along the three years. Among those environmental predictors, uncorrelated ones were selected in order to avoid any identification confusion in Maxent due to their redundancies. Correlations between all possible pairs of predictors were calculated across the background grid range using Spearman's determination coefficient (Fig. S3, Supplementary Material), together with Principal Components Analysis (Fig. S4, Supplementary Material). When the determination coefficient of a pair of predictors exceeded 0.7, and their first and second components on the PCA plot were close, one of the two predictors was discarded (Sillero et al. 2021, Fig. S4). Among the 13 potential predictors, 10 were not correlated: five geomorphological predictors (Depth, Slope, Rugosity, Eastness and Northness) and five hydrological predictors (Umax, annual average SST, Tmin, Tmax and annual average Turbidity). For each of these 10 predictors, their statistical distribution was assessed across the five species' presence location sets, across the entire background grid and on restrictions of the background grid to each of the eight sites to evidence potential sampling bias. All data processing was performed using in-house © Matlab routines.

### Diagnostics of model fitting

Model best fit was the average of five replicated probability maps fitted using all predictors under the cross-validate option. For each replicated fit, the presence data was partitioned randomly (random seed option) into a training dataset used to fit the model (75 % of presence data) and a testing dataset (25 % of presence data). AUC on the training and testing datasets and the relative contribution of each predictor to the replicated fit were calculated to assess the quality of the fit. The contribution of each predictor in defining each of the five replicated fits was ranked using a jackknife cross-validation procedure (Phillips et al. 2006). To do so, for each of the five presence data partitions, replicated fits were constructed using either a predictor on its own or all predictors except the one in question. Gains and losses of AUC compared to the

AUC with all predictors were calculated with and without each predictor, respectively.

Another ranking of the contribution of each predictor to a model was performed using a permutation procedure. For each predictor in turn, the values of that predictor in the training presence and background data were randomly permuted. The model was then re-evaluated on the permuted data, and the resulting drop in training AUC indicated that the model depended heavily on that predictor. Values were normalized to give percentages.

Finally, predictor response curves showed how the predicted probability of presence varied with the environmental predictor holding all other environmental predictors at their mean sample value.

### Sensitivity of model fitting to potential sampling bias when reducing presence number and their geographical extent

Sampling bias can appear when some environmental conditions are preferentially sampled without respecting the proportion of their availability (Merow et al. 2013). Sampling bias can be the result of geographical or temporal preference in presence exploration or limited presence exploration. Such bias is evident for rare species whose niche is limited in space and requires intensive exploration to be observed. For more widely distributed species, sampling bias is less evident. Among the five gorgonian species, we selected the two species with highest occurrence: *E. singularis* (255 presence data) and *L. sarmentosa* (62 presence data), and assessed the sensitivity of the AUC when decreasing the number of presence data in the training dataset (while still keeping the full presence dataset for testing). To do so, we followed two procedures. In the first procedure, the model was fitted using three training datasets, gradually reducing the number of presence data to an evenly distributed selection of 50 %, 20 % and 10 % of the initial training dataset. In the second procedure, the model was fitted using a reduced number of presence data in the training datasets, corresponding to geographical restrictions to six sites for *E. singularis* (CVM, AGD, ARES, AGM, PMCB and PNPC) and two sites for *L. sarmentosa* (CVM and AGD). For both species, a training dataset including only presence recorded within MPAs was tested as well. Finally, sampling bias was also tested by testing if adding the geographical coordinates as predictors altered the model fitting.

### Model performance assessment

The performance of two models, one using the full presence dataset and the other only using the presence observed in MPAs, was compared. Model performance was assessed by contrasting the habitat suitability predictions with observed presences and absences for each of the five species (the 444 grid cells subset mentioned above).

For each species, regional habitat suitability maps were built after binarizing the mean of the five replicated probability maps. This was done by applying a threshold probability above which the habitat was considered as suitable and below which it was considered as unsuitable. Due to varying prevalence among the five species, the threshold probability
was not set to a fixed value but instead was determined using one of the 11 possible thresholding procedures advised in Maxent (Phillips *et al.* 2017). Deciding which threshold to use is a major consideration when minimizing false positives and false negatives. Its choice, however, is often arbitrary in studies (Liu *et al.* 2005). Among these threshold procedures, the 10percentile \_training\_ threshold (maximum error rate of 10 % on the presence points of the learning data) was arbitrarily selected because it displayed the highest threshold values when compared to other thresholding procedures, resulting in the less extended and hence more conservative habitat suitability maps (Fig. S5, Supplementary Material).

# RESULTS

## Species presence observations and sampling bias assessment

In the study area, the observed spatial distribution varied among the five species. Some species were distributed in a limited number of clustered locations: *C. rubrum* with 39 presences in only two sites (CVM and PMBC, Fig. 2), *P. clavata* with 61 presences in three sites (CVM, PMBC and PNPC, Fig. 2) and *E. cavolini* with 87 presences distributed in three sites (PNPC, PMCB and ARES, Fig. 2). *Leptogorgia sarmentosa* was detected everywhere except in PNPC and AGM (Fig. 2). Finally, *E. singularis* was the most widespread species with 479 presences, distributed in all eight sites (Fig. 2). The species most frequently associated together were *E. singularis* and *L. sarmentosa* (17.9 %). Boxplots that represented the environmental predictors range, showed marked differences among the eight sites forming the background seascape (Fig. 3). Depth range was doubled in four out of the eight sites, extending down to 40 m depth (CVM, AGD, PMCB and PNPC, Fig. 3D). Larger slope (Fig. 3E) and rugosity (Fig. 3F) were found in two (PNPC and CVM) out of the eight sites, while seabed orientation (eastness and northness) displayed a similar range of variability in all eight sites (Fig. 3G and H). Higher annual average SST (Fig. 3B) characterized one of the eight sites (PNPC) while Tmin at the seabed (Fig. 3J) clearly separated two sites (PMCB and PNPC) where temperature at the seabed did not fall below 11 °C from the other six sites where minimum temperature was below 10 °C. Lower annual average turbidity (Fig. 3C) characterised three out of the eight sites (CVM, PMCB and PNPC). Finally, larger Umax at the

seabed was frequently found in five out of the eight sites (LEU, VLR, AGD, ARES, AGM, Fig. 3A).

Boxplots that represented environmental predictors range at species presence locations, indicated some differences between species as well. Umax at the seabed was consistently higher at presence locations of *C. rubrum* than of any other species (Fig. 3A). Annual average SST was consistently lower at presence locations of *C. rubrum*, *E. singularis* and *L. sarmentosa* than of *P. clavata* and *E. cavolini* (Fig. 3B). Annual average turbidity was higher at presence locations of *E. singularis* and *L. sarmentosa* (Fig. 3C), while slope was slightly higher at presence locations of *P. clavata* than at other species presence locations (Fig. 3E).

## Environmental predictors contribution to habitat suitability prediction

Different environmental predictors contributed to the Maxent model best fit of each gorgonian species in the Gulf of Lion (when only using the presence locations found in the background grid, Fig. 2). For *C. rubrum*, annual average turbidity and Umax at the seabed contributed more than 90 % to the suitable habitat prediction (both in percent contribution and permutation percentage, Table I). Response curves for these two predictors indicated species preference for current speeds greater than 0.25 m.s⁻¹ and annual average turbidity lower than 449 m-1 (Fig. S7, Supplementary Material), suggesting that the ecological niche of *C. rubrum* extends beyond the variation range of those predictors within the background grid of the study. Annual average turbidity and Umax were the predictor that, when absent, most decreased the AUC in the model, and Umax at the seabed contributed the most to the AUC in stand-alone mode, followed by depth and rugosity (Jackknife's test, Fig. S6, Supplementary Material). For *P. clavata*, Tmin at the seabed, depth, and Umax combined contributed more than 70 % of the prediction of suitable habitat (percentage contribution, Table I). However, Tmin at the seabed contributed only 6.3 % in permutation percentage (Table I), reflecting the similarity in the range of values of Tmin in the background and in
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 2.- Presence location (Blue cross) and regional habitat suitability maps predicted using Maxent with all predictors except geographical ones (Latitude / Longitude) (Grid red). The habitat suitability map was calculated by applying a threshold value to obtain binary maps (0-1) from Cloglog output. A: close up PNPC, B: close up LEU, C: close up CVM, D: close up AGD, E: close up AGM and ARES, F: close up VLR, G: close up PMCB. Maps for : (1) *C. rubrum*, (2) *P. clavata*, (3) *E. cavolini*, (4) *E. singularis*, (5) *L. sarmentosa*.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 2a: Presence location (Blue cross) and regional habitat suitability maps predicted using Maxent for *E. cavolini* with all predictors including geographical ones (Latitude / Longitude) (Grid red). The habitat suitability map was calculated by applying a threshold to obtain binary maps (0-1) from Cloglog output. A: close up PNPC, B: close up LEU, C: close up CVM, D: close up AGD, E: close up AGM and ARES, F: close up VLR, G: close up PMCB. Maps for : (1) *C. rubrum*, (2) *P. clavata*, (3) *E. cavolini*, (4) *E. singularis*, (5) *L. sarmentosa*.</figcaption></figure>
<table><caption>Table I.- Significant environmental predictors identified for each species by Maxent without geographical predictors (Latitude / Longitude). Percent contribution and permutation importance</caption><thead><tr><th></th><th>Variables</th><th>Percent contribution (%)</th><th>Permutation importance (%)</th></tr></thead><tbody><tr><td rowspan="2">C. rubrum</td><td>Umax</td><td>70.5</td><td>43.3</td></tr><tr><td>Annual average turbidity</td><td>22.9</td><td>49.8</td></tr><tr><td rowspan="4">P. clavata</td><td>Tmin</td><td>31.1</td><td>6.3</td></tr><tr><td>Depth</td><td>26.1</td><td>35.7</td></tr><tr><td>Umax</td><td>14.6</td><td>18.8</td></tr><tr><td>Annual average turbidity</td><td>6.2</td><td>19.5</td></tr><tr><td rowspan="3">E. cavolini</td><td>Tmin</td><td>32.5</td><td>45.2</td></tr><tr><td>Rugosity</td><td>30.3</td><td>33</td></tr><tr><td>Depth</td><td>20.3</td><td>13.4</td></tr><tr><td rowspan="4">E. singularis</td><td>Rugosity</td><td>25.9</td><td>13.6</td></tr><tr><td>Annual average SST</td><td>22.3</td><td>24.6</td></tr><tr><td>Depth</td><td>14.6</td><td>26.9</td></tr><tr><td>Annual average turbidity</td><td>13.7</td><td>12.8</td></tr><tr><td rowspan="6">L. sarmentosa</td><td>Depth</td><td>37</td><td>6</td></tr><tr><td>Rugosity</td><td>21.6</td><td>20.2</td></tr><tr><td>Annual average SST</td><td>10.7</td><td>10.3</td></tr><tr><td>Annual average turbidity</td><td>10.7</td><td>17.8</td></tr><tr><td>Tmin</td><td>5.2</td><td>18.9</td></tr><tr><td>Tmax</td><td>4.2</td><td>20.2</td></tr></tbody></table>

*P. clavata* presence sites (Fig. 3B). Conversely, annual average turbidity was ranked fourth in the percentage of contribution but second in percentage of permutation (19.5 %, Table I), reflecting the narrow range of annual average turbidity in *P. clavata* presence sites compared to the annual average turbidity range in the background (Fig. 3C).

Response curves of the four predictors indicated species preference for Tmin larger than 8.2 °C, water depth greater than 24 m, current greater than 0.14 m.s-1, and annual average turbidity less than 357 m-1 (Fig. S8, Supplementary Material). All response curves displayed a plateau when they reached high values, suggesting the
species can distribute beyond 50 m water depth and into stronger currents. Depth was the predictor that most decreased the AUC when it was not included in the model (Jackknife's test, Fig. S6).

For *E. singularis*, rugosity, annual average SST, depth contributed similarly to the suitable habitat prediction (both in percent contribution and permutation percentage, Table I), summing to 75 %. Response curves indicated that *E. singularis* had a preference for annual average SST ranging from 15.7 and 17.5 °C, water depth ranging from 4 to 40 m, with a characteristic niche-like optimal range of values and rugosity larger than 1 (Fig. S10, Supplementary Material). Annual average SST was the predictor that most decreased the AUC when removed from the model and contributed the most to the AUC in stand-alone mode (Jackknife's test, Fig. S6).

For *L. sarmentosa*, depth, rugosity, annual average SST and turbidity contributed the most to the suitable habitat prediction, summing up to 80 % in percent contribution. However, in permutation percentage, Tmin and Tmax at the seabed replaced annual average SST and depth, whereas together with rugosity and annual average, turbidity contributed to 77 % (Table I). Response curves indicated that a preferential habitat for *L. sarmentosa* would be in areas with water depth greater than 15 m, a temperature below 22.8 °C, in terrain with rugosity between 1 and 5 and with an annual average turbidity of between 204 and 600 m-1 (Fig. S11, Supplementary Material). No predictor had a significant influence on AUC when they were absent from the model. Annual average SST, depth and Tmax contributed the most to stand-alone model AUC (Jackknife's test, Fig. S6).

In summary, annual average turbidity was identified as a significant environmental predictor for four species. Maximum flow speed at the seabed was the primary predictor for two species out of five. Among geomorphological predictors, rugosity and depth contributed to suitable habitat modelling, while no orientation predictors did. Despite rugosity being a significant predictor for all species, ru-

gosity response curves did not indicate much discrimination between the species, as any rugosity larger than 1 was favourable. Minimum flow speed at the seabed did not contribute to the suitable habitat modelling of any of the five gorgonian species.

### Model performance and assessment of the predictions

Average AUC varied from 0.805 to 0.985 (less than 5 % standard deviation) across the five species, indicating good to very good fit when learning and testing against replicated random partitions of the presence dataset (Table II). Maps of suitable habitats that were obtained after applying a thresholding occurrence probability to the five replicated models were highly congruent, with less than 10 % of locations where the five replicated model predictions disagreed for all species except *E. singularis* (20 %). Maps of suitable habitats derived from the model best fit demonstrated correct predictions when compared to both presence and absence observations for all species (Fig. 2). However, model accuracy was lower for *L. sarmentosa* and *E. singularis* (70 % and 73 %, respectively) compared to *P. clavata* and *C. rubrum* (92 % for both, Table II). In all species, TSS differed significantly from a random contingency table (Chi-square, P < 0.0001), with values larger than 0.8 for *C. rubrum* and *P. clavata* and larger than 0.4 for *L. sarmentosa* and *E. cavolini* (Table II). For *E. singularis*, TSS was only 0.33 due to low specificity (43 % of absence correctly predicted, Table II), reflecting the large extension of the suitable habitat predicted for the species.

### Assessing the influence of the number and location of presence observations used for model training on model fit goodness

Sensitivity of model AUC to the number and location of presence observations used for model training was carried out on *E. singularis* and *L. sarmentosa* because these species have the widest spatial coverage and the greatest number of presences. For *E. singularis* and *L. sarmentosa*, model AUC remained within the range of variation of the five replicated model AUCs when the training dataset was reduced by 90 % and 80 %, respectively,
A

<figure><img src="image_9.png" /></figure>

C

<figure><img src="image_10.png" /></figure>

E

<figure><img src="image_11.png" /></figure>

B

<figure><img src="image_12.png" /></figure>

D

<figure><img src="image_13.png" /></figure>

F

<figure><img src="image_14.png" /></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 3.- Boxplot of the 10 non-correlated predictors (except latitude and longitude) in the 8 sites (CVM, LEU, VLR, AGD, ARES, AGM, PMCB, PNPC), at species presence locations and in the 100 m x 100 m background grid (A: Umax, B: Annual average SST, C: Annual average turbidity, D: Depth, E: Slope, F: Rugosity, G: Northness, H: Eastness, I: Tmax, J: Tmin).</figcaption></figure>

as long as the presence observations of the training dataset were spread in all eight sites (Fig. 4). In contrast, model AUC decreased sharply when the number of presence observations was reduced by the same amount but restricted to a single site only (Fig. 4). When the number of presence observations was restricted to the highly protected zones of the area, model AUC decreased to 0.7 for *E. singularis* and 0.85 for *L. sarmentosa*. This decrease was, however, less than if the same number of presences had been taken from a single site, although slightly more than if it was evenly distributed in all eight sites. When the number of presence observations was restricted to the highly protected zones of the area, model AUC decreased only for *P. clavata* and

*L. sarmentosa* (0.93 and 0.88 respectively, Table II) while accuracy decreased for *P. clavata* and *E. singularis* (0.89 and 0.59 respectively, Table II). Adding latitude and longitude as environmental predictors, the mean test AUC was unchanged for *E. singularis*, *P. clavata* and *C. rubrum*, increased for *E. cavolini* (0.955) and decreased for *L. sarmentosa* (0.890). For *L. sarmentosa*, adding geographical predictors altered predicted distributions, resulting in the loss of a site where the species has been observed (ARES, Fig. S12, Supplementary Material). In contrast, the suitable habitat map for *E. cavolini* reduced considerably. Suitable habitat area was removed in CVM and reduced in AGD and AGM with improvement of accuracy.
<table><caption>Table II.- Assessment of the performance of model and predicted suitability maps trained on 80 % of presence data with all predictors except geographical ones (Lat / Lon) : AUC test average and Standard deviation of 5 models, Accuracy, Sensitivity: true positif, Specificity: true absence. Values in Italic between parentheses corresponds to model trained on sampling in AMPs only.</caption><thead><tr><th></th><th>AUC test</th><th>Standard<br>deviation of<br>AUC test</th><th>Accuracy</th><th>Sensitivity</th><th>Specificity</th><th>TSS</th></tr></thead><tbody><tr><td><i>C. rubrum</i></td><td>0.985(0.989)</td><td>0.005(0.019)</td><td>0.92(0.94)</td><td>0.90(0.65)</td><td>0.92(0.95)</td><td>0.82(0.60)</td></tr><tr><td><i>P. clavata</i></td><td>0.974(0.936)</td><td>0.007(0.023)</td><td>0.92(0.89)</td><td>0.88(0.88)</td><td>0.92(0.89)</td><td>0.81(0.77)</td></tr><tr><td><i>E. cavolini</i></td><td>0.932(0.965)</td><td>0.045(0.037)</td><td>0.82(0.94)</td><td>0.90(0.10)</td><td>0.81(1)</td><td>0.71(0.1)</td></tr><tr><td><i>E. singularis</i></td><td>0.805(0.932)</td><td>0.034(0.034)</td><td>0.73(0.59)</td><td>0.91(0.44)</td><td>0.43(0.85)</td><td>0.33(0.29)</td></tr><tr><td><i>L. sarmentosa</i></td><td>0.911(0.886)</td><td>0.027(0.073)</td><td>0.70(0.78)</td><td>0.77(0.57)</td><td>0.68(0.83)</td><td>0.45(0.40)</td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>Fig. 4.- Assessment of the AUC for *Eunicella singularis* (white triangle) and *Leptogorgia sarmentosa* (dark cercle) varying the number of presence data (down to 7 % of the initial training dataset) and their location in the training dataset while keeping the full presence data in the testing dataset.</figcaption></figure>

cy and TSS (0.96 and 0.86 respectively, Fig. 2a)

## DISCUSSION

The observations between 10 and 50 m depth across eight rocky sites spanning 450 km along the Gulf of Lion (NW Mediterranean Sea) confirmed that between 2013 and 2020 the hierarchy in the occurrence of the five gorgonian species was similar to the one previously reported

in the same area (Weinberg 1979a, Linares et al. 2008, Gori et al. 2011a). *Eunicella singularis* was the most abundant gorgonian species with a wide continuous spatial distribution. *Leptogorgia sarmentosa* was the second most abundant gorgonian species, with scattered distribution. *Paramuricea clavata* and *Corallium rubrum* displayed highly patchy wwdistributions and their presence was observed in three rocky sites only. *Eunicella cavoli*
ni was the least abundant species, with marginal presence in a few spots in three rocky sites. Occasional associations between the five gorgonian species were observed but were not frequent, similar to previously reported between 0 and 70 m deep south of the Gulf of Lion (Gori *et al.* 2011a).

Depth is a primary factor that has been used as a basis for the classification of marine benthic communities strata in the shallow environment as both light availability and agitation at sea decrease with depth (Perez 1961, Odum 1971). In the present study, depth was identified as a significant predictor for all species distribution except *C. rubrum*, but with different responses according to the species. For *P. clavata*, model predictions corroborated a clear preference for deeper habitats as already mentioned in the literature (Rossi 1959, Boavida *et al.* 2016, Fourt *et al.* 2017). Similarly, model predictions for *E. singularis* which is the only one of the five species hosting *Symbiodiniacea* (Forcioli *et al.* 2011) indicated a preference for shallow areas which might be related to more abundant light (between 4 and 32 m). Nevertheless, a different morphotype of the same species with extremely low densities of *Symbiodiniacea* has been reported beyond 50 m depth (Theodor 1969, Gori *et al.* 2011b). Moreover, limited gene flow between the shallow and deep morphotypes was reported at local scale (Costantini *et al.* 2016). These observations suggest that care should be taken when delineating habitat suitability for this species, as different morphotypes may likely display different niches, hence related to different environmental predictors. For *L. sarmentosa*, depth response excluded depth less than 15 meters while the species has been observed in shallow harbour areas (Betti *et al.* 2018). For *E. cavolini*, depth response only excluded very shallow environments (shallower than 5 m) in the present study area while in the locality of Kornati (Adriatic Sea), an *E. cavolini* population was found at 5 m depth (Sini *et al.* 2015). For *C. rubrum*, depth was not identified as a significant predictor, despite the fact that the species preference for low light has been largely reported (Laborel 1960, Weinberg 1979b, Bianconi *et al.* 1988). The presence of *C. rubrum*

at shallow depth might, instead, be related to the smaller size of this species, compared to the other four gorgonian species, which would allow the colonization of small dark crevices (Laborel 1961, Rossi *et al.* 2008). In any case, despite the fact that depth is easier to measure than light, direct measurement or modelling of light intensity should be undertaken to improve the delineation of the species niche (Betti *et al.* 2018).

Seabed orientation is another proxy for light availability that has been used in several studies at small spatial scales (Glasby 2000, Knott *et al.* 2004). However, in the present study, orientation was never indicated as a significant predictor. This was, most probably, due to the fact that orientation and light availability are not correlated over a large spatial scale.

Beyond depth, the five gorgonian species have been associated with different preferences in geomorphological features. *C. rubrum* is mainly found in overhangs or in caves (Cau *et al.* 2016). *P. clavata* and *E. cavolini* are reported to dwell in steep rocky coralligenous bottoms, such as drop-offs, and they are often found in association (Carugati *et al.* 2021). *L. sarmentosa* have been associated to the detritic-muddy benthic community, characterized by boulders and abundance of organic matter (Carpine & Grasshoff 1975, Mistri & Ceccherelli 1993). *E. singularis* has been mainly reported on well-lit horizontal and sub-horizontal rocky bottoms (Gori *et al.* 2011a, Weinberg & Weinberg 1979). The latter two species have also been reported on gently sloping rocky bottoms (Gori *et al.* 2011a). High-resolution terrain mapping thanks to airborne and satellite remote sensing has been an essential data source for determining species habitats when using predictive habitat modeling in terrestrial environments (Wilson *et al.* 2007). However, high-resolution seabed mapping over large areas is still rare. Vessel-mounted multibeam acoustic sounders can provide seabed maps with a vertical resolution at the cm scale and horizontal resolution at the tens of cm scale, but only over a small surface (Dolan *et al.* 2008). The cost of acquiring such data has limited its implementation, which is why predictive ha-
bitat modelling studies at sea have focused on
specific environments such as canyons (Guinan
et al. 2009, Bargain et al. 2018, Lo Iacono et al.
2018). The present study aimed to fill a gap in
predictive habitat modelling for shallow marine
environments by leveraging advances in bathymet-
try remote sensing using airborne laser imaging
between 2011 and 2015 (LITTO3D, www.diffu-
sion.shom.fr), which enabled mapping shallow
areas along hundreds of kms of coastline, with
submetric spatial resolution similar to that of v
essel-mounted acoustic sounders. Yet, our choice
to average geomorphological predictors over a 20
m x 20 m square slab to match the decametric scale
at which the species presence was detected, no lon-
ger enabled to detect rapid terrain variations with
a single variability predictor. This is a frequent
dilemma when quantifying geomorphological va-
riations (Wilson et al. 2007), and it explains why
we kept two variability predictors (slope and rugo-
sity), as in combination they have the potential to
differentiate flat tilted bottom (low slope and large
rugosity) from trenches (large slope and low ru-
gosity) and walls (large slope and large rugosity).
However, model predictions did not indicate slope
as a significant environmental predictor for any of
the five species, while steep slope is often assu-
med to favor benthic invertebrates colonisation, as
the sediment deposition on flat surfaces could limit
their recruitment (Laborel & Vacelet 1961, Bian-
coni et al. 1988, Glasby 2000). In particular, walls,
which have often been described as the ecological
niche of *P. clavata*, may not actually be (Cocito et
al. 2002, Boavida et al. 2016, Ponti et al. 2019).
This result is consistent with the presence of very
big colonies of *P. clavata* on horizontal rocks in the
deep environment (Grinyó et al. 2016). Reversely,
rugosity differentiated the suitable habitat of *L.
sarmentosa*, which is mainly found at a moderate
rugosity, and that of other two species (*E. singula-
ris* and *E. cavolini*) which can be found on a wider
range of rugosity. Yet, the LIDAR resolution did
not enable the detection of crevices smaller than
0.1 m in which *C. rubrum* colonies are often found
(True 1970, Virgilio et al. 2006). This may explain
why rugosity was not identified as a significant
predictor for this latter species in the present study.

In addition to geomorphological pre-
dictors, hydrological predictors have been put
forward as factors regulating the presence of
gorgonians. Among these, turbidity is one of the
trickiest ones. Sea water turbidity remotely sensed
from satellites can indicate sedimentation due to
river inputs or to resuspension by waves (Liew
et al. 2009, Gohin et al. 2020). High sedimenta-
tion is considered a limiting factor for most gor-
gonians, potentially due to the burial and/or clo-
gging of recently settled polyps, except for L.
sarmentosa (Otero et al. 2017). The preference of
L. sarmentosa for turbid seawater would be related
to a morphological adaptation, with higher polyp
density per branch and finer branches, which re-
duced surface would limit sediment accumulation
in turbid environmental conditions, compared
to the other four gorgonian species (Cocito et al.
2002, Rossi et al. 2011). In the present study, an-
nual average turbidity was identified as a signifi-
cant predictor for all species, except E. cavolini. As
expected, C. rubrum and P. clavata were limited
to a lower range of turbidity than L. sarmentosa,
which displayed a well delimited range of interme-
diate turbidity preference (Weinberg 1979a, Mis-
tri 1995). For E. singularis, the response to turbi-
dity was unclear with tolerance to both high and
low annual average turbidity, which may reflect
the presence of two morphotypes in the 10-50 m
depth range (with high or low densities of Sym-
biodiniacea, Gori et al. 2012). Indeed, preference
for turbid environments can be expected in hete-
rotrophic suspension feeders and might be related
to food requirements while preference for less tur-
bid environments can be expected in autotrophic
species relying on photosynthesis. Apart from the
broad differentiation between autotrophy and he-
terotrophy dominance in E. singularis, available
environmental predictors did not enable distingui-
shing between living (zooplankton) and non-living
 detritic) food preferences of the other four hete-
rotrophic gorgonian species. Anyhow, P. clavata,
L. sarmentosa and C. rubrum display only
subtle differences in the proportions of
their trophic sources (Ribes et al. 1999,
Ribes et al. 2003, Tsounis et al. 2006).
Sedimentation can be locally reduced by large
flow speed values (Hiscock 1983), another factor suggested as structuring the distribution of gorgonians (Leversee 1976, Weinberg 1979a) and deep-water corals (White et al. 2005). The five studied gorgonian species, being passive suspension feeders, are expected to have a preference for moderate to strong flow speed values (Gori et al. 2011a, La Riviere et al. 2021). The present study identified large flow speed as a significant predictor for *C. rubrum* and *P. clavata* only. For *P. clavata*, maximum flow speed at the seabed was a more significant predictor than slope and roughness to describe the spatial distribution of the species, pointing to the need of adding hydrodynamical predictors to avoid confusion of factors when analysing the drivers of species distribution. Indeed, flow speed at the seabed is not commonly used as a predictor in ecological niche models (but see Bargain et al. 2018). It should be more systematically tested, taking care of using hydrodynamical simulations at adequate spatial resolution near the coast (Sciascia et al. 2022). Finally, temperature appeared as a significant environmental predictor for four of the five species, but never highlighting a limitation due to large maximum temperature at the seabed, while temperature positive anomalies during heat waves have been related to mass mortalities in *E. singularis* (Ezzat et al. 2013), *E. cavolini* (Turicchia et al. 2018) and *P. clavata* (Coma et al. 2006, Bally & Garrabou 2007, Garrabou et al. 2009). Although *L. sarmentosa* is found even inside ports, its tolerance to high water temperature (up to 29°C) was not evidenced either (Betti et al. 2018, personal observation). Actually, temperature predictors for *E. singularis* and *L. sarmentosa* (annual average SST < 17.5°C) and *P. clavata* and *E. cavolini* (Tmin > 8°C) reflect large-scale geographical separation between the colder Gulf of Lion sites and the warmer PNPC and PMCB sites. In summary, identifying the environmental predictors for marine species still faces limitations and will require both increasing the number of presence observations to deeper areas and including more explicitly factors such as local temperature, light, turbidity, sedimentation, chlorophyll-a and flow speed. To this aim, the use of high-resolution flow

and primary production simulations would be an invaluable tool for ecological studies, as they allow projecting the distribution of the future suitable habitats under climate change scenarios and test the deep refugia hypothesis (Bongaerts et al. 2010).

Nevertheless, despite the above-mentioned limitations, Maxent performance was reasonable for *E. singularis* and very good for the other four species (according to AUC rating, Hanley & McNeil, 1982, Sainz-Villegas et al. 2022). Accuracy of predicted suitable habitat maps was good to excellent for all species (according to TSS rating, Araujo et al. 2011). The good sensitivity of the models for *E. singularis* and *L. sarmentosa* with few but regularly distributed observations of presence, confirms, on the one hand, that spatial autocorrelation was avoided and that large AUC values could not be attributed to model overfitting (Radosavljevic & Anderson 2013) and, on the other hand, that a low number of presence data is sufficient to obtain a good performance of the models as long as a wide geographical area is encompassed, to avoid local specificities in the training process. Moreover, in shallow environments, characterised by a diversity of environmental conditions and strong human impacts, we advocate for the use of highly protected zones (such as the Gulf of Lion MPAs) as sentinel sites, to allow for the monitoring of the factors shaping the benthic species distribution in coastal areas. Indeed, gorgonians like other benthic species are known to be sensitive to mechanical impacts such as mooring (Broad et al. 2020), bottom fishing (Bavestrello et al. 1997, Betti et al. 2020), scuba diving (Coma & Zabala 1994) and even direct exploitation (i.e., *C. rubrum* for jewellery, Tsounis et al. 2010, 2013). The small number of presence observations obtained from MPA species inventories, that were repeated over time, provided enough data to understand the habitat preferences of these species. In the present study, true absence and presence data used to assess model accuracy were defined based on repeated inventories obtained between 2013 and 2020 and avoided bias due to transient population accidents (*E. singularis* was observed in 2013 but was absent in 2019 in some locations). Given the
slow population dynamics of gorgonians, impac-
ted populations may take several years to recover
(Linares et al. 2007). The reiteration of sampling
over a decade is thus essential, in order to assess
the ENMs performance. In addition, an observa-
tion strategy that is focused on a limited number of
locations would improve the measurement of pre-
dictors by equipping sentinel sites with adequate
sensors (e.g., seawater temperature measurement
networks such as T-MEDNet, Garrabou et al. 2019).

For all species, the suitable habitat predic-
tion area was greater than the original inventory
range. As explained above, suitable habitat is, by
definition, expected to extend beyond actual pre-
sence observations and repeated observations in
time are key to the fine tuning of true absence data
and model accuracy. However, in the present stu-
dy, *E. cavolini* displayed very reduced presence
distributions compared with its predicted suitable
habitat. Moreover, geographical predictors impro-
ved the prediction. A phylogenetic study suggested
that differentiation between *E. singularis* and *E.
cavolini* may have resulted from isolation in diffe-
rent refugia during the last glaciation (Aurelle *et*
al. 2017). In such a hypothesis, the absence of *E.
cavolini* from the Gulf of Lion (except one locality
in ARES) could be attributed to a quicker re-colo-
nisation of the shallow suitable habitat by *E. singu-
laris* after the last glaciation, pointing to relatively
recent geographical vicariance. Since then, the two
species seem to have evolved towards diverging en-
vironmental niches as their predictors differ. Howev-
ver, beyond the definition of a suitable habitat based
on environmental predictors (potential niche), ac-
tual species distribution (realized niche) ultimately
depends on efficient colonisation or demographic
connectivity, which is a complex process that ope-
rates at multiple spatial and temporal scales, resul-
ting from species demography (population size,
structure, fecundity) and larval dispersal (larval
release timing, ocean currents, pelagic larval du-
ration, larval behaviour). Scenarios using spatially
explicit metapopulation modelling can help tes-
ting the respective role of the suitable habitat and
demographic connectivity in explaining observed
species spatial distribution (Guizien et al. 2014).

To conclude, habitat suitability models

differentiated the niche of *P. clavata*, *C. rubrum*, *E. singularis*, *E. cavolini* and *L. sarmentosa* using seven predictors (maximum flow speed at the bottom, annual average turbidity, depth, minimum and maximum temperature at the bottom, sea surface temperature and rugosity). The spatial distribution of the suitable habitat was predicted with excellent accuracy and sensitivity for the five species in the Gulf of Lion. These models can be used to complement the delineation of *P. clavata*, *C. rubrum*, *E. singularis*, *E. cavolini* and *L. sarmentosa* suitable habitat in the shallow rocky bottoms in other regions along the Mediterranean basin and guide marine spatial planning, including the establishment of new MPAs. However, observations in the mesophotic and aphotic environments should be added to the models developed here for the shallow rocky bottoms and should be updated to incorporate more direct predictors (i.e., light, sedimentation rate, organic matter content) in order to avoid the use of proxies such as depth and turbidity.

**ACKNOWLEDGEMENTS.**- This work was funded by the French National Program LITEAU IV of the Ministère de l'Ecologie et de l'Environnement Durable under project RocConnect - Connectivité des habitats rocheux fragmentés du Golfe du Lion (PI, K. Guizien, Project Number 12-MUTS-LITEAU-1-CDS-013), by the EC Interreg MArittimo program under project IMPACT - IMpatto Portuale su aree marine protette: Azioni Cooperative Transfrontaliere (PI, M. Magaldi, Ecology Task Leader, K. Guizien, CUP B12F17000370005), and by the Agence de l'Éau Rhône-Mé diterranée-Corse under project ICONE - Impacts actuels et potentiels de la CONnectivité Ecologique ajoutée par les récifs artificiels sur la biodiversité fixée des substrats durs du Golfe du Lion (PI, K. Guizien, AAP 2016). The authors gratefully acknowledge the helpful assistance during sam- pling of R. Bricout, B. Hesse, L. Lescure, J. wwC. Roca, and the staff of the Aire Marine Protégée Agathoise, the Réserve Naturelle Nationale Marine de Cerbère-Banyuls and the Parc Marin de la Côte Bleue. Thanks to the SHOM (https://dif- fusion.shom.fr/) for the LITTO3D lidar data and the OFB, the MPAs mentionned above for the multibeam sonar data and the habitat maps. Thanks to NASA Goddard Space Flight Center, Ocean Ecology Laboratory, Ocean Biology Proces- sing Group, (2018): Sea-viewing Wide Field-of-view Sen- sor (SeaWiFS) Ocean Color Data, NASA OB.DAAC. doi: 10.5067/ORBVIEW-2/SEAWIFS/L2/OC/2018. Accessed on 06/06/2017 and 20/05/2022 for AQUAMODIS data.

Authors contributions

SB and KG conceived the study, SB and LB carried out sam-
pling and SB did the statistical analysis. All contributed to
manuscript writing.

**Conflict of interest disclosure**

The authors of this preprint declare that they have no financial conflict of interest with the content of this article.

**REFERENCES**

Achu AL, Gopinath G, Surendran U 2021. Landslide susceptibility modelling using deep-learning and machine-learning methods-A study from southern Western Ghats. *In* International India Geoscience and Remote Sensing Symposium (InGARSS), Ahmedabad, India, 6-10 December 2021: 360-364. https://doi.org/10.1109/InGARSS51564.2021.9792034

Allouche O, Tsoar A, Kadmon R 2006. Assessing the accuracy of species distribution models: Prevalence, kappa and the true skill statistic (TSS): Assessing the accuracy of distribution models. *J Appl Ecol* 43 (6): 1223-1232. https://doi.org/10.1111/j.1365-2664.2006.01214.x

Araújo MB, Alagador D, Cabeza M, Nogués-Bravo D, Thuiller W 2011. Climate change threatens European conservation areas: Climate change threatens conservation areas. *Ecol Lett* 14: 484-492. https://doi.org/10.1111/j.1461-0248.2011.01610.x

Aurelle D, Pivotto ID, Malfant M, Topcu NE, Masmoudi MB *et al.* 2017. Fuzzy species limits in Mediterranean gorgonians (Cnidaria, Octocorralla): Inferences on speciation processes. *Zool Scr* 46: 767-778. https://doi.org/10.1111/zsc.12245

Bailey M 2010. Ecosystem-based Management for the Oceans, 44 (2): 304-305. *In* Island Press. McLeod K, & Heather L, (EDS), Washington, DC, USA. https://doi.org/10.1017/S0030605310000244

Bally M, Garrabou J 2007. Thermodependent bacterial pathogens and mass mortalities in temperate benthic communities: A new case of emerging disease linked to climate change. *Glob Change Biol* 13 (10): 2078-2088. https://doi.org/10.1111/j.1365-2486.2007.01423.x

Bargain A, Foglini F, Pairaud I, Bonaldo D, Carniel S *et al.* 2018. Predictive habitat modeling in two Mediterranean canyons including hydrodynamic variables. *Prog Oceanogr* 169: 151-168. https://doi.org/10.1016/j.pocean.2018.02.015

Bavestrello G, Cerrano C, Zanzi D, Cataneo-Vietti R 1997. Damage by fishing activities to the Gorgonian coral *Paramuricea clavata* in the Ligurian Sea. *Conserv Mar Freshw Ecosyst* 7 (3): 253-262. https://doi.org/10.1002/(SICI)1099-0755(199709)7:3<253::AID-AQC243>3.0.CO;2-1

Betti F, Bavestrello G, Bo M, Ravaneti G, Enrichetti F *et al.* 2020. Evidences of fishing impact on the coastal gorgonian forests inside the Portofino MPA (NW Mediterranean Sea). *Ocean Coast Manag* 187: 105. https://doi.org/10.1016/j.ocecoaman.2020.105105

Betti F, Bo M, Bava S, Faimali M, Bavestrello G 2018. Shallow-water sea fans: The exceptional assemblage of *Leptogorgia sarmen-tosa* (Anthozoa: Gorgoniidae) in the Genoa harbour (Ligurian Sea). *Eur Zool J* 85 (1): 290-298. https://doi.org/10.1080/24750263.2018.1494219

Bianconi CH, Rivoire G, Stiller A, Boudouresque CF 1988. Le corail rouge *Corallium rubrum* (Lamarck) dans la réserve naturelle de Scandola. *Trav scient Parc Nat Reg réser nat Corse* 16: 1-83.

Boavida J, Assis J, Silva I, Serrão EA 2016. Overlooked habitat of a vulnerable gorgonian revealed in the Mediterranean and Eastern Atlantic by ecological niche modelling. *Sci Rep* 6 (1): 36460. https://doi.org/10.1038/srep36460

Bongaerts P, Ridgway T, Sampayo EM, Hoegh-Guldberg O 2010. Assessing the 'deep reef refugia' hypothesis: Focus on Caribbean reefs. *Coral Reefs* 29 (2): 309-327. https://doi.org/10.1007/s00338-009-0581-x

Bramanti L, Vielmini I, Rossi S, Stolfa S, Santangelo G 2011. Involvement of recreational scuba divers in emblematic species monitoring: The case of Mediterranean red coral (Corallium rubrum). *J Nat Conserv* 19 (5): 312-318. https://doi.org/10.1016/j.jnc.2011.05.004

Breiman L 2001. Random Forests. *Machine Learning* 45: 5-32. https://doi.org/10.1023/A:1010933404324

Briton F, Cortese D, Duhaut T, Guizien K 2018. High-resolution modelling of ocean circulation can reveal retention spots important for biodiversity conservation. *Aquat Conserv Mar Freshw Ecosyst* 28 (4): 882-893. https://doi.org/10.1002/aqc.2901
Broad A, Rees MJ, Davis AR 2020. Anchor and chain scour as disturbance agents in benthic environments: Trends in the literature and charting a course to more sustainable boating and shipping. *Mar Pollut Bull* 161: 111683. https://doi.org/10.1016/j.marpolbul.2020.111683

Canessa M, Bavestrello G, Trainito E 2023. *Leptogorgia sarmentosa* (Anthozoa: Octocoralia) in NE Sardinia (Mediterranean Sea): distribution and growth patterns. *Mar Biodivers* 53: 13. https://doi.org/10.1007/s12526-022-01313-0

Carpine C, Grasshoff M 1975. Les Gorgonaires de la Méditerranée. *Bull Inst Océano Monaco*. Vol 71(1430): 1-14.

Carrasco PC 2010. Nugget effect, artificial or natural? *J South Afr Inst Min Metall* 110: 299-305.

Carugati L, Moccia D, Bramanti L, Cannas R, Follesa MC et al. 2022. Deep-Dwelling Populations of Mediterranean *Corallium rubrum* and *Eunicella cavolini*: Distribution, Demography, and Co-Occurrence. *Biology* 11 (2): 333. https://doi.org/10.3390/biology11020333

Carugati L, Moccia D, Cau A, Bramanti L, Cannas R et al. 2021. *Corallium rubrum* and *Eunicella cavolini*: Distribution, population structure and co-occurrence in the deep Mediterranean Sea. *In* International Workshop on Metrology for the Sea, Learning to Measure Sea Health Parameters (MetroSea), Reggio Calabria, Italy, 04-06 October 2021: 95-99. https://doi.org/10.1109/MetroSea52177.2021.9611596

Cau A, Bramanti L, Cannas R, Follesa MC, Angiolillo M et al. 2016. Habitat constraints and self-thinning shape Mediterranean red coral deep population structure: Implications for conservation practice. *Sci Rep* 6 (1):23322. https://doi.org/10.1038/srep23322

Çinar ME, Gönulal O, Özturk B 2018. Wanted dead or alive: *Corallium rubrum* (Cnidaria: Anthozoa) on the coasts of Turkey. *Cah biol mar* 59: 175-179. https://doi.org/10.21411/CBM.A.1043F293

Cocito S, Bedulli D, Sgorbini S 2002. Distribution patterns of the sublittoral epibenthic assemblages on a rocky shoal in the Ligurian Sea (NW Mediterranean). *Sci Mar* 66 (2): 175-181. https://doi.org/10.3989/scimar.2002.66n2175

Cohen J 1960. A coefficient of agreement

of nominal scales. *Educat Psychol Meas* 20: 37-46.

Coma R, Linares C, Ribes M, Diaz D, Garrabou J et al. 2006. Consequences of a mass mortality in populations of *Eunicella singulAIRIS* (Cnidaria: Octocorallia) in Menorca (NW Mediterranean). *Mar Ecol Prog Ser* 327:51-60. doi.org/10.3354/meps327051

Coma R, Zabala M 1994. Efecto de la frequentacion sobre las poblaciones de *ParamuRICEA clavata* en las illes Medes. *In* VIII Simposio Iberico de Estudios del Bentos Marino, Blanes, Spain. Universitat de Barcelona: 170-171.

Costantini F, Gori A, Lopez-González P, Bramanti L, Rossi S et al. 2016. Limited genetic connectivity between gorgonian morphotypes along a depth gradient. *PLoS ONE* 11 (8): e0160678.

doi:10.1371/journal.pone.0160678

Cúrdia J, Monteiro P, Afonso CML, Santos MN, Cunha MR, Gonçalves JMS 2013. Spatial and depth-associated distribution patterns of shallow gorgonians in the Algarve coast (Portugal, NE Atlantic). *Helgol Mar Res* 67: 521-534. doi.org/10.1007/s10152-012-0340-1

Di Camillo CG, Ponti M, Bavestrello G, Krzelj M, Cerrano C 2018. Building a baseline for habitat-forming corals by a multi-source approach, including Web Ecological Knowledge. *Biodivers Conserv* 27 (5): 1257-1276. doi.org/10.1007/s10531-017-1492-8

Dolan MFJ, Grehan AJ, Guinan JC, Brown C 2008. Modelling the local distribution of cold-water corals in relation to bathymetric variables: Adding spatial context to deep-sea video data. *Deep Sea Res Part Oceanogr Res Pap* 55 (11): 1564-1579. doi.org/10.1016/j.dsr.2008.06.010

Elith J, Graham C, Anderson R, Dudík M, Ferrier S et al. 2006. Novel methods improve prediction of species' distributions from occurrence data. *Ecography* 29 (2): 129-151. https://doi.org/10.1111/j.2006.0906-7590.04596.x

Elith J, Phillips SJ, Hastie T, Dudík M et al. 2011. A statistical explanation of MaxEnt for ecologists. *Divers Distrib* 17 (1): 43-57. doi.org/10.1111/j.1472-4642.2010.00725.x

Evans IS 1980. An integrated system of terrain analysis and slope mapping. *In*:
tgr Syst Terrain Anal Slope Mapp 36: 274-295.

Ezzat L, Merle PL, Furla P, Buttler A, Ferrier-Pagès C 2013. The Response of the Mediterranean Gorgonian *Eunicella singularis* to Thermal Stress Is Independent of Its Nutritional Regime. *PLoS ONE* 8 (5): e64370. https://doi.org/10.1371/journal.pone.0064370

Forcioli D, Merle PL, Caligara C, Ciosi M, Muti C et al. 2011. Symbiont diversity is not involved in depth acclimation in the Mediterranean sea whip *Eunicella singularis*. *Mar Ecol Prog Ser* 439: 57-71. https://doi.org/10.3354/meps09314

Fournier M, Goujard A, Pérez T, Chevaldonné P 2017. Guide de la faune profonde de la mer Méditerranée : Explorations des roches et canyons sous-marins des côtes françaises. Muséum national d'Histoire naturelle, Paris. *Patrimoines naturels*, Tome 75: 184 p.

Garrabou J, Coma R, Bensoussan N, Baly M, Chevaldonné P et al. 2009. Mass mortality in Northwestern Mediterranean rocky benthic communities: Effects of the 2003 heat wave. *Glob Change Biol* 15 (5): 1090-1103. https://doi.org/10.1111/j.1365-2486.2008.01823.x

Gerovasileiou V, Smith CJ, Sevastou K, Papadopoulou N, Dailianis T et al. 2019. Habitat mapping in the European Seas - is it fit for purpose in the marine restoration agenda? *Mar Policy* 106: Article 103521, 10.1016/j.marpol.2019.103521.

Gili JM, Coma R 1998. Benthic suspension feeders: Their paramount role in littoral marine food webs. *Trends Ecol Evol* 13: 316-321. https://doi.org/10.1016/S0169-5347(98)01365-2

Glasby TM 2000. Surface composition and orientation interact to affect subtidal epibiota. *J Exp Mar Biol Ecol* 248 (2): 177-190. https://doi.org/10.1016/S0022-0981(00)00169-6

Gohin F, Bryère P, Lefebvre A, Sauriau PG, Savoye N et al. 2020. Satellite and In Situ Monitoring of Chl-a, Turbidity, and Total Suspended Matter in Coastal Waters: Experience of the Year 2017 along the French Coasts. *J Mar Sci Eng* 8 (9): 665. https://doi.org/10.3390/jmse8090665

Gori A, Bramanti L, López-González P, Thoma, JN, Gili JM et al. 2012. Characterization of the zooxanthellate and azooxanthellate morphotypes of the Mediterranean gorgonian

*Eunicella singularis*. *Mar Biol* 159 (7): 1485-1496. https://doi.org/10.1007/s00227-012-1928-3

Gori A, Rossi S, Berganzo E, Paretus JL, Dale MRT et al. 2011a. Spatial distribution patterns of the gorgonians *Eunicella singularis*, *Paramuricea clavata*, and *Leptogorgia sarmentosa* (Cape of Creus, Northwestern Mediterranean Sea). *Mar Biol* 158 (1): 143-158. https://doi.org/10.1007/s00227-010-1548-8

Gori A, Rossi S, Linares C, Berganzo E, Orejas C et al. 2011b. Size and spatial structure in deep versus shallow populations of the Mediterranean gorgonian *Eunicella singularis* (Cap de Creus, northwestern Mediterranean Sea). *Mar Biol* 158: 1721-1732.

López-González PJ 1993. Taxonomia y zoogeografia de los antozoos del Estrecho de Gibraltar y areas proximas.Ph.D.thesis,UniversidaSevilla,Spain.

Grasshoff M 1977. Die gorgonarien des oestlichen nordatlantik und des mittelmeeres. Die familie Paramuriceidae (Cnidaria, Anthozoa). *METEOR Forschungsergebnisse, herausgegeben von der deutschen Forschungsgemeinschaft* 27: 5-76.

Grinyó J Gori, A, Ambroso S, Purroy A, Calatayud C et al. 2016. Diversity, distribution and population size structure of deep Mediterranean gorgonian assemblages (Menorca Channel, Western Mediterranean Sea). *Prog Oceanogr* 145: 42-56. https://doi.org/10.1016/j-pocean.2016.05.001

Guinan J, Brown C, Dolan MFJ, Grehan AJ 2009. Ecological niche modelling of the distribution of cold-water coral habitat using underwater remote sensing data. *Ecol Inform* 4 (2): 83-92. https://doi.org/10.1016/j.ecoinf.2009.01.004

Guisan A, Thuiller W, Zimmermann NE 2017. Habitat suitability and distribution models. With Applications in R (Ecology, Biodiversity and Conservation). In Cambridge: Cambridge University Press. doi: 10.1017/9781139028271

Guizien K, Belharet M, Moritz C, Guarnini JM 2014. Vulnerability of marine benthic metapopulations: implications of spatially structured connectivity for conservation practice in the Gulf of Lions (NW Mediterranean Sea). *Divers Distrib*
20: 1392-1402. https://doi.org/10.1111/ddi.12254

Guizien K, Bramanti L, Roca JC, Hesse B, Blouet S et al. 2022. Database of sea fans spatial distribution in the northwestern Mediterranean Sea. https://www.seanoe.org/data/00750/86176/ https://doi.org/10.17882/86176

Hanley JA, McNeil BJ 1982. The meaning and use of the area under a receiver operating characteristic (ROC) curve. *Radiology* 143 (1): 29-36. https://doi.org/10.1148/radiology.143.1.7063747

Hernandez PA, Franke I, Herzog SK, Pachec V, Paniagua L et al. 2008. Predicting species distributions in poorly-studied landscapes. *Biodivers Conserv* 17 (6): 1353-1366. https://doi.org/10.1007/s10531-007-9314-z

Hernandez PA, Graham CH, Master LL, Albert DL 2006. The effect of sample size and species characteristics on performance of different species distribution modeling methods. *Ecography* 29 (5): 773-785. https://doi.org/10.1111/j.0906-7590.2006.04700.x

Hiscock K 1983. Water movement. The ecology of the shallow sublittoral benthos. *In* Sublittoral ecology, Earl R & Erwin DG (Eds), Oxford: Clarendon Press : 58-96.

Horta e Costa B, Claudet J, Franco G, Erzini K, Caro A et al. 2016. A regulation-based classification system for Marine Protected Areas (MPAs). *Mar Policy* 72, 192-198. https://doi.org/10.1016/j.marpol.2016.06.021

Jarnevich CS, Holcombe TR, Thomas CC, Frid L, Olsson A 2015. Simulating long-term effectiveness and efficiency of management scenarios for an invasive grass. *AIMS Environ Sci* 2 (2): 427-447. https://doi.org/10.3934/environsci.2015.2.427

Knittweis L, Aguilar R, Alvarez H, Borg JA, Evans J 2016. New depth record oh the red coral *Corallium rubrum* for the Mediterranean. *Rapp Com Intern Explor Scie Mer Méd* 41: 467.

Knott NA, Underwood AJ, Chapman MG, Glasby TM 2004. Epibiota on vertical and on horizontal surfaces on natural reefs and on artificial structures. *J Mar Biol Assoc UK* 84 (6): 1117-1130. https://doi.org/10.1017/S0025315404010550h

La Rivière M, Michez N, Delavenne J, Andres S, Fréjefond C et al. 2021. Fiches descripives des biocénoses benthiques de Méditerranée.

UMS PatriNat (OFB-CNRS-MNHN): 660 p.

Laborel J 1960. Contribution à l'étude directe des peuplements benthiques scaphiles sur substrats rocheux en Méditerranée. *Rec trav stat Mar Endoume* 33 (20): 117-173.

Laborel J 1961. Le concrétionnement algal «coralligène» et son importance géomorphologique en Méditerranée. *Rec trav stat mar Endoume* 23: 37-60.

Labore J, Vacelet J 1961. Répartition bio-nomique du *Corallium rubrum* (Lmck) dans les grottes sous-marines. *Rapp réun commiss intern explor scient Méd Monaco* 16: 464-469.

Legendre P, Dale MRT, Fortin MJ, Gurevitch J, Hohn M et al. 2002. The consequences of spatial structure for the design and analysis of ecological field surveys. *Ecography* 25 (5): 601-615. https://doi.org/10.1034/j.1600-0587.2002.250508.x

Leversee GJ 1976. Flow and feeding in fan-shaped colonies of the gorgonian coral, *Leptogorgia. Biol Bull* 151 (2): 344-356. https://doi.org/10.2307/1540667

Liew SC, Saengtuksin B, Kwoh LK 2009. Monitoring turbidity and suspended sediment concentration of coastal and inland waters using satellite data. *In* International Geoscience and Remote Sensing Symposium, IGARSS 2009, July 12-17, 2009, University of Cape Town, Cape Town, South Africa: 837-839. https://doi.org/10.1109/IGARSS.2009.5418225

Linares C, Coma R, Garrabou J, Díaz D, Zabala M 2008. Size distribution, density and disturbance in two Mediterranean gorgonians: *Paramuricea clavata* and *Eunicella sin-gularis. J Appl Ecol* 45 (2): 688-699. https://doi.org/10.1111/j.1365-2664.2007.01419.x

Linares C, Doak, DF, Coma R, Díaz D, Zabala M 2007. Life history and viability of a long-lived marine invertebrate: the oc-tocoral *Paramuricea clavata*. *Ecology* 88 (4):918-928. https://doi.org/10.1890/05-1931

Liu C, Berry, PM, Dawson TP, Pearson RG 2005. Selecting thresholds of occurrence in the prediction of species distributions. *Ecography* 28 (3): 385-393. https://doi.org/10.1111/j.0906-7590.2005.03957.x

Lo Iacono C, Robert K, Gonzalez-Villanueva.
va R, Gori A, Gili JM et al. 2018. Predicting cold-water coral distribution in the Cap de Creus Canyon (NW Mediterranean): Implications for marine conservation planning. *Prog Oceanogr* 169: 169-180. https://doi.org/10.1016/j-pocean.2018.02.012

McCullagh P 1983. Generalized linear models (2nd ed). Routledge, New York: 532 p. doi: 10.1201/9780203753736

McFadden CS, Van Ofwegen LP, Quatrini AM 2022. Revisionary systematics of Octocorallia (Cnidaria: Anthozoa) guided by phylogenomics. *Bull Soc Syst Biol* 1(3): 79 p. https://doi.org/10.18061/bssb.v1i3.8735

Merckx B, Steyaert M, Vanreusel A, Vincx M, Vanaverbeke J 2011. Null models reveal preferential sampling, spatial autocorrelation and overfitting in habitat suitability modelling. *Ecol Model* 222 (3): 588-597.

Merow C, Smith, MJ, Silander Jr JA 2013. A practical guide to Maxent for modeling species' distributions: What it does, and why inputs and settings matter. *Ecography* 36 (10): 1058-1069. https://doi.org/10.1111/j.1600-0587.2013.07872.x

Mistri M 1995. Population Structure and Secondary Production of the Mediterranean Octocoral *Lophogorgia ceratophyta* (L, 1758). *Mar Ecol* 16 (3), 181-188. https://doi.org/10.1111/j.1439-0485.1995.tb00404.x

Mistri M, Ceccherelli VU 1993. Growth of the Mediterranean Gorgonian *Lophogorgia ceratophyta* (L, 1758). *Mar Ecol* 14 (4), 329-340. https://doi.org/10.1111/j.1439-0485.1993.tb00004.x

Nikuradse J, 1950. Stromungsgesetze in rauhen Rohren. Forschungsheft 361: 65 p.

Odum EP 1971. An understanding of ecological succession provides a basis for resolving man's conflict with nature. *Science* Vol 164, 262-270.

Otero MM, Numa C, Bo M, Orejas C, Garrabou Jet al. 2017. Overview of the conservation status of Mediterranean Anthozoa. International Union for Conservation of Nature (IUCN), Malaga, Spain: 73 p. https://doi.org/10.2305/IUCN.CH.2017.RA.2.en

Pearson RG 2007. Species' Distribution Modeling for Conservation Educators and Practitioners. Synthesis America Museum of Natural History. American Museum of Natural History 50 p.

Perez JM 1961. Océanographie bio-

logique et biologie marine. Tome premier: la vie benthique. Les Presses Universitaires de France, Collection Euclide, Paris, 541 p.

Phillips SJ, Anderson RP, Schapire RE 2006. Maximum entropy modeling of species geographic distributions. *Ecolological Modelling* 190 (3-4): 231-259. https://doi.org/10.1016/j.ecolmodel.2005.03.026

Phillips SJ, Anderson R, Dudík M, Schapire R, Blair M 2017. Opening the black box: An open-source release of Maxent. *Ecography* 40: 887-893. https://doi.org/10.1111/ecog.03049

Phillips SJ, Dudík M 2008. Modeling of species distributions with Maxent: New extensions and a comprehensive evaluation. *Ecography* 31: 161-175.

Ponti M, Turicchia E, Costantini F, Gori A, Bramanti L et al. 2019. Mediterranean gorgonian forests: Distribution patterns and ecological roles. In *Proceedings of the 3rd Mediterranean symposium on the conservation of coralligenous & other calcareous Bio-concretions, Antalya, Turkey, 15-16 January 2019*. SPA/RAC, Tunis: 4-14.

Radosavljevic A, Anderson R 2013. Making better Maxent models of species distributions: Complexity, overfitting and evaluation. *J Biogeoogr* 41: 629-643. https://doi.org/10.1111/jbi.12227

Ribes M, Coma R, Gili J 1999. Heterogeneous feeding in benthic suspension feeders: the natural diet and grazing rate of the temperate gorgonian *Paramuricea clavata* (Cnidaria: Octocorallia) over a year cycle. *Mar Ecol Prog Ser* 183: 125-137. https://doi.org/10.3354/meps183125

Ribes M, Coma R, Rossi S 2003. Natural feeding of the temperate asymbiotic octocoral-gorgonian *Leptogorgia sarmentosa* (Cnidaria: Octocorallia). *Mar Ecol Prog Ser* 254: 141-150. https://doi.org/10.3354/meps254141

Rossi L 1959. Le specie di Eunicella (Gorgonaria) del golfe di Genova. *Ann Mus ci vico storia nat Giacomo Doria* 71: 203-225.

Rossi S, Gili JM, Garrofé X 2011. Net negative growth detected in a population of *Leptogorgia sarmentosa*: Quantifying the biomass loss in a benthic soft bottom-gravel gorgonian. *Mar Biol* 158 (7): 1631-1643. https://doi.org/10.1007/s00227-011-1675-x

Rossi S, Tsounis G, Orejas C, Padrón T, Gili,
JM, Bramanti L *et al*. 2008. Survey of deep-dwelling red coral (*Corallium rubrum*) populations at Cap de Creus (NW Mediterranean). *Mar Biol* 154 (3), 533-545. https://doi.org/10.1007/s00227-008-0947-6

Sainz-Villegas S, De La Hoz CF, Juanes JA, Puente A 2022. Predicting non-native seaweeds global distributions: The importance of tuning individual algorithms in ensembles to obtain biologically meaningful results. *Front Mar Sci* 9: 1009808.

Salomidi M, Smith C, Katsanevakis S, Papayotidis P, Papathanassiou V 2009. Some observations on the structure and distribution of gorgonian assemblages in the Eastern Mediterranean Sea: 242-245. *In Actes du 1er Symposium sur la Conservation du Coralligène et autres bioconcrétions de Méditerranée, Tabarka, 16-17 Janvier 2009.*

Sciascia R, Guizien K, Magaldi MG 2022. Larval dispersal simulations and connectivity predictions for Mediterranean gorgonian species: sensitivity to flow representation and biological traits. *ICES J Mar Sci* 79: 2043-2054. https://doi.org/10.1093/icesjms/fsac135

Sillero N, Arenas-Castro S, Enriquez-Urzelai U, Vale CG, Sousa-Guedes D *et al.* 2021. Want to model a species niche? A step-by-step guideline on correlative ecological niche modelling. *Ecol Model* 456: 109671. https://doi.org/10.1016/j.ecolmodel.2021.109671

Sini M. Kipson S, Linares C, Koutsoubas D, Garrabou J 2015. The Yellow Gorgonian *Eunicella cavolini*: Demography and Disturbance Levels across the Mediterranean Sea. *PLoS ONE* 10(5): e0126253: 19 p.

Theodor J 1969. Contribution à l'étude des gorgones (Viii): *Eunicella stricta* aphyta sous espèce nouvelle sans zooxanthelles, proche d'une espéce normalement infestée par ces algues. Documents faunistiques et écologiques. *Vie Milieu*: 635-637.

Topçu EN, Öztürk B 2015. Composition and abundance of octocorals in the Sea of Marmara, where the Mediterranean meets the Black Sea. *Sci Mar* 79: 125-135. https://doi.org/10.3989/scimar.04120.09A

True MA 1970. Étude quantitative de quatre peuplements sciaphiles sur substrat rocheux dans la région marsellaise. *Bull Inst Océano Monaco* 69 (1401): 1-48.

Tsounis G, Rossi S, Bramanti L, Santangelo G 2013. Management hurdles in sustainable harvesting of *Corallium rubrum*. *Mar Policy* 39: 361-364. https://doi.org/10.1016/j.marpol.2012.12.010i

Tsounis G, Rossi S, Gili JM *et al*. 2006. Population structure of an exploited benthic cnidarian: the case study of red coral (*Coralium rubrum* L.). *Mar Biol* 149: 1059-1070. https://doi.org/10.1007/s00227-006-0302-8

Tsounis G, Rossi S, Grigg R, Santangelo G, Bramanti L *et al*. 2010. The exploitation and conservation of precious corals. *Ocean. Mar Biol Annu Rev* 48: 161-212.

Turicchia E, Abbiati M, Sweet M, Ponti M 2018. Mass mortality hits gorgonian forests at Montecristo Island. *Dis Aquat Organ* 131(1):79-85. https://doi.org/10.3354/dao03284

Vafidis D 2009. First record of *Leptogorgia sarmentosa* (Octocorallia: Gorgoniidae) from the eastern Mediterranean Sea. *Mar Biodivers Rec* 2: e17. Doi:10.1017/S1755267208000195

Virgilio M, Airoldi L, Abbiati M 2006. Spatial and temporal variations of assemblages in a Mediterranean coralligenous reef and relationships with surface orientation. *Coral Reefs* 25 (2): 265-272. https://doi.org/10.1007/s00338-006-0100-2

Vissenaekens E, Guizien K, Durrieu de Madron X, Pairaud I, Leredde Y, Puig P, Bourrin F 2023. Accuracy of high resolution coastal flow speed simulations during and outside of wind, wave and stratification events (Gulf of Lion, NW Mediterranean). *J Mar Syst* 239 (103845): 1-13 https://doi.org/10.1016/j.jmarsys.2022.103845

Weinberg S 1978. The minimal area problem in invertebrate communities of Mediterranean rocky substrata. *Mar Biol* 49 (1), 33-40. https://doi.org/10.1007/BF00390728

Weinberg S 1979a. Autecology of Shallow-Water *Octocorallia* from Mediterranean Rocky Substrata, I. The Banyuls Area. *Bijdr Tot Dierkd* 49(1), 1-15. https://doi.org/10.1163/26660644-04901001

Weinberg S 1979b. The Light-Dependent Behaviour of Planula Larvae of *Eunicella singularis* and *Corallium rubrum* and its Implication for Octocorallian Ecology. *Bijdr Tot Dierkd* 49 (1): 16-30. https://doi.org/10.1163/26660644-04901002

Weinberg S, Weinberg F 1979. The Life
Cycle of a Gorgonian: *Eunicella Singularis* (Es-
per, 1794). *Contrib Zool Bijdr Tot Dierkd* 48: 127-
140. https://doi.org/10.1163/26660644-04802003

White M, Mohn C, De Stigter H, Mot-
tram G 2005. Deep-water coral development as
a function of hydrodynamics and surface pro-
ductivity around the submarine banks of the
Rockall Trough, NE Atlantic. *In* Cold-Water
Corals and Ecosystems. A Freiwald & J M Ro-
berts (Eds.) Springer, Berlin, Heidelberg: 503-
514. https://doi.org/10.1007/3-540-27673-4_25

Wilson MFJ, O'Connell B, Brown C, Gui-
nan JC, Grehan AJ 2007. Multiscale terrain analy-
sis of multibeam bathymetry data for habitat map-
ping on the continental slope. *Mar Geod* 30 (1-2): 3-35. https://doi.org/10.1080/01490410701295962

Wisz MS, Hijmans RJ, Li J, Peterson AT, Graham CH, Guisan A 2008. Effects of sample size on the performance of species distribution models. *Divers Distrib* 14: 763-773.

Zibrowius H, Marques VM, Grasshoff M 1984. La répartition du *Corallium rubrum* dans l'atlantique (Cnidaria: Anthozoa: Gorgonaria). *Téthys* 11(2): 163-170.

Zweig MH, Campbell G 1993. Receiver-operating characteristic (ROC) plots: a fundamental evaluation tool in clinical medicine. *Clin Chem* 39 (8): 561-57.

Received February 9, 2024

Accepted March 19, 2024

Associate editor: F. Lartaud
# SUPPLEMENTARY MATERIALS

## ECOLOGICAL NICHE MODELLING OF FIVE GORGONIAN SPECIES WITHIN THE SHALLOW ROCKY HABITAT OF THE FRENCH MEDITERRANEAN COAST

S. BLOUETI,<fnref n="2" />*, L. BRAMANTI<fnref n="1" />, K. GUIZIEN<fnref n="1" />

## MATERIALS AND METHODS

### Background grid, features and regularization

The background grid is a set of georeferenced locations where probability of species presence will be predicted (a prior with uniform species prevalence is defined) using known environmental predictors values. Probability distribution of environmental predictors in the background grid is compared to probability distribution of environmental predictors in presence locations in order to build response curves. Hence, Maxent background grid actually defines the range of environmental conditions on which features will be fitted in order to maximize the species entropy probability distribution, resulting in different features and Maxent predictions according to background extent (Phillips et al. 2006).

Various feature classes can be combined in Maxent to build response curves from environmental predictors. Theoretically, a niche model is expected to be unimodal, pointing at species optimum, which will be best represented by quadratic features. Sometimes the niche range, however, might be truncated in the presence observations. In this situation, linear or hinge features will be more appropriate (Merow et al. 2013). Multiplying feature class possibilities may lead to confusion because complex features created by Maxent are often highly correlated. From a statistical perspective, when Maxent is used to infer environmental drivers of species distribution, feature classes should be selected with parsimony. From a machine-learning perspective, testing as many feature classes as possible will increase the predictive accuracy of the presences. In this study, we adopted an in-between perspective, selecting the default «auto-features» mode of Maxent. This mode increases the tested feature classes to the number of presence observations with all feature classes for more than 80 presence observations (Phillips & Dudík 2008). This mode enabled Maxent to accommodate the different number of presence observations, which varied from 16 for *C. rubrum* to 255 for *E. singularis*. Regularisation dampens predictors coefficients when predictors variance is high at the presence locations. Regularisation multiplier was varied from 0.1 to 5 and the AUC of model fit varied by less than 1 %. Hence, it was set to its default value of 1.

### Output type

The ROR is the relative probability that a cell is contained in a collection of presence samples, assuming individuals have been randomly sampled in the environmental space. A logistic transformation of ROR enables one to predict the probability of presence-only, assuming the probability of presence in the suitable habitat, namely the prevalence, is 50 % (Phillips & Dudík 2008). Due to the difficulty of checking this assumption, the safer interpretation of Maxent predictions was using ROR as a habitat suitability index (Merow et al. 2013). As an alternative, the cloglog value calculated as $1-\exp(-c \cdot r)$ has been proposed in order to estimate the probability of presence (c is the exponential of the entropy of the Maxent distribution and r is the ROR value, Phillips et al. 2006). In the present study, however, model fit was not sensitive to the species prevalence value in the background grid. The latter varied from 20 to 80 % and resulted in average AUC (Area Under the Curve) variation of 0.6 %. Therefore, prevalence was kept to its default value of 50 % and the cloglog output type was selected.

<footnote n="1">¹ CNRS-Sorbonne Université, Laboratoire d'Ecogéochimie des Environnements Benthiques, LECOB, Observatoire Océanologique de Banyuls Sur Mer, 66650 Banyuls sur Mer, France.</footnote>

<footnote n="2">² Ville d'Agde, Aire marine protégée de la côte agathoise, 34300 Agde, France.</footnote>

<footnote>* Corresponding author: sylvain.blouet@ville-agde.fr</footnote>
## Definitions of the seven geomorphological predictors

All predictors were computed on the 20 m x 20 m slab located at the center of the 100 m x 100 m grid cell.

* Local bottom depth was the average of bathymetric data within the slab.

* A local digital terrain model (DTM) was defined by fitting all bathymetric

al data within the slab with a bivariate quadratic equation (1) (Evans 1980):

$$Z = aX^2 + bY^2 + c XY + dX + e Y + f \quad (1)$$

with Z being the bottom depth deviation from the mean bottom depth and X, Y the cartesian distance to the center of the background grid cell and a, b, c, d, e and f are the fitting parameters that defines the local DTM and vary at each background grid cell. These parameters are subsequently used to calculate the slope, roughness, rugosity, ruggedness and orientation. Slope was defined as the average slope of the local DTM:

$$S = \arctan(\sqrt{(d^2+e^2)})$$

* Roughness is length scale (in meter) defined as the difference between the maximum and minimum depth in the slab which is consistent with the hydraulic roughness parameter used to estimate seabed friction by Nikuradse (1933).

* Rugosity is a topological metric (no unit) defined as the ratio between the surface of the local DTM and the slab surface. By definition, rugosity is larger or equal to 1.

* Terrain ruggedness index was defined as the average of absolute deviation of the local DTM with respect to the mean depth.

* Eastness and northness orientation were defined as the sinus and cosinus of the DTM aspect, the latter being defined as:

$$Aspect = atan(d/e)$$

<figure><img src="image_17.png" /><figcaption>Fig. S1.- Sampling location and bathymetry in each of the eight sites. A: close up PNPC, B: close up LEU, C: close up CVM, D: close up AGD, E: close up AGM and ARES, F: closeup VLR, G: close up PMCB.</figcaption></figure>
<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>
Fig. S2.- Empirical semivariograms describing the spatial variability of each of the five gorgonian species in the four sites (ARES, AGD, CVM, PMCB) where the species was present. Semivariograms were computed using population spatial density measured in the 696 evenly spaced geo-referenced locations covering the hard bottom habitat and clustering them by distance classes ranging from 500 meters to 10 kms, with a step of 500 meters increment. The nuggwwet of the semivariogram was computed using replicated measurements made at less than 20 meters distance, around each geo-referenced location (Guizien et al. 2022).

<figure><img src="image_23.png" /></figure>

Fig. S3.- Spearman determination coefficient among the 13 initials environmental predictors including latitude and longitude.
<figure><img src="image_24.png" /><figcaption>Fig. S4.- Principal Components Analysis on the 13 initials environmental predictors including latitude and longitude.</figcaption></figure>
A

<figure><img src="image_25.png" /></figure>

Threshold

B

<figure><img src="image_26.png" /></figure>

Threshold
C

E.cavolinii

Fixed_cumulative_value_1
Fixed_cumulative_value_5
Minimum-trainingpresence

Fixed_cumulative_value_10
Minimum-training presenc
Minimum training presence
Minimum training plus specific
MinimumCumulativeValue_5
MinimumCumulativeValue_10
MinimumCumulativeValue_5
MinimumCumulativeValue_10

Fixed_cumulative_value_10
Minimumtrainingpresence
Minimumtraining plus specific
Equal-training sensitivity and spec
Maximum_test_sensitivity plus specific
Maximum_test_sensitivity_and_spec
Maximum_test_sensitivity plus spec
Maximum_test_sensitivity_and_specification

Fixed_cumulative_value_10
Minimumtrainingpresence
Minimumtraining plus specificity
Equal-training sensitivity and specificity
Maximum_test_sensitivity plus specificity
Maximum_test_sensitivity_and_spec
Maximum_test_sensitivity plus specific
Maximum_test_sensitivity_and_specification

Balance-training omission
Balance training omission
Equate_entropy_of_thresholded_andof thresholded and origina
Equate_entropy_of_thresholded_andof thresholded and origina

Threshold

D

E. singularis

Fixed_cumulative_value_1
Fixed_cumulative_value_5
Minimumtrainingpresence

Fixed_cumulative_value_10
Minimumtrainingpresence
Minimumtraining plus specific
Equal-training sensitivity and specific
Maximum_test_sensitivity plus specific
Maximum_test_sensitivity_and_spec

Fixed_cumulative_value_10
Minimumtrainingpresence
Minimumtraining plus specificity
Equal-training sensitivity and specificity
Maximum_test_sensitivity plus specificity
Maximum_test_sensitivity_and_spec

Fixed_cumulative_value_10
Minimumtrainingpresence
Minimumtraining plus specificity
Equaltraining sensitivity and specific
Maximum_test_sensitivity plus specificity
Maximum_test_sensitivity_and_specification

Balance-training omission
Balance training omission
Equate_entropy_of_thresholded_and
Equate_entropy_of_thresholded_and_Original_distributions

Threshold
<figure><img src="image_27.png" /><figcaption>Fig. S5.- Threshold proposed by Maxent to predict suitability maps for various indices and for each species. Modelling with all predictors except geographical ones (Latitude / Longitude). A: *C. rubrum*, B: *P. clavata*, C: *E. cavolini*, D: *E. singularis*, E: *L. sarmentosa*.</figcaption></figure>
<figure><img src="image_28.png" /></figure>

A

<figure><img src="image_29.png" /></figure>

B

<figure><img src="image_30.png" /></figure>

C
<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>

Fig S6.- Maxent AUC applying two Jackknife procedures across all predictors except geographical ones (Latitude / Longitude). A: *C. rubrum*, B: *P. clavata*, C: *E. cavolini*, D: *E. singularis*, E: *L. sarmentosa*. Without variable (light blue), with only variable (blue), with all variables (red). Variable is the same as predictor.
<figure><img src="image_33.png" /><figcaption>Fig. S7.- Response curves for the two significant predictors in the Maxent model for habitat suitability of *C. rubrum*: A: Umax, B: Annual average turbidity. The red lines indicate the mean values, while blue areas the standard deviation limits, resulting from 5 cross-validation model runs in cloglog output.</figcaption></figure>

<figure><img src="image_34.png" /></figure>
<figure><img src="image_35.png" /><figcaption>Fig. S8.- Response curves for the four significant predictors in the Maxent model for habitat suitability of *P. clavata*. A: Tmin, B: Depth, C: Umax, D: Annual average turbidity. The red lines indicate the mean values, while blue areas the standard</figcaption></figure>
<figure><img src="image_36.png" /></figure>

Fig. S9: Response curves for the three significant predictors in the Maxent model for habitat suitability of *E. cavolini*. A: Tmin, B: Rugosity, C: Depth. The red lines indicate the mean values, while blue areas the standard deviation limits, resulting from 5 cross-validation model runs in cloglog output.
<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /></figure>

<figure><img src="image_39.png" /></figure>
<figure><img src="image_40.png" /></figure>

Fig. S10.- Response curves for the four significant predictors in the Maxent model for habitat suitability of *E.singularis*. A: Rugosity, B: Annual average SST, C: Depth, D: Annual average turbidity. The red lines indicate the mean values, while blue areas the standard deviation limits, resulting from 5 cross-validation model runs in cloglog output.

<figure><img src="image_41.png" /></figure>
<figure><img src="image_42.png" /></figure>

<figure><img src="image_43.png" /></figure>

<figure><img src="image_44.png" /></figure>
<figure><img src="image_45.png" /><figcaption>Fig. S11.- Response curves for the six significant predictors in the Maxent model for habitat suitability of *L. sarmentosa*. A: Rugosity, B: Annual average SST, C: Annual average turbidity, D: Tmin, E: Maximum bottom temperature Tmax, F: Depth. The red lines indicate the mean values, while blue areas the standard deviation limits, resulting from 5 cross-validation model runs in cloglog output.</figcaption></figure>

<table><caption>Table S1.- Significant environmental predictors identified for each species by Maxent with geographical predictors (Latitude / Longitude). Percent contribution and permutation importance.</caption><thead><tr><th></th><th>AUC</th><th>Variables</th><th>Percent contribution (%)</th><th>Permutation importance (%)</th></tr></thead><tbody><tr><td rowspan="4">C. rubrum</td><td rowspan="4">0.9628</td><td>Umax</td><td>56.5</td><td>21.7</td></tr><tr><td>Latitude</td><td>26.4</td><td>0</td></tr><tr><td>Annual average turbidity</td><td>6.9</td><td>34.8</td></tr><tr><td>Tmin</td><td>4.5</td><td>32.6</td></tr><tr><td rowspan="4">P. clavata</td><td rowspan="4">0.9731</td><td>Tmin</td><td>30.5</td><td>14.6</td></tr><tr><td>Depth</td><td>27.9</td><td>32.2</td></tr><tr><td>Umax</td><td>12.8</td><td>13</td></tr><tr><td>Annual average turbidity</td><td>3.4</td><td>14.8</td></tr><tr><td rowspan="3">E. cavolini</td><td rowspan="3">0.9491</td><td>Rugosity</td><td>27.5</td><td>52.9</td></tr><tr><td>Depth</td><td>20.5</td><td>14.8</td></tr><tr><td>Latitude</td><td>5.8</td><td>18</td></tr><tr><td rowspan="3">E. singularis</td><td rowspan="3">0.8133</td><td>Latitude</td><td>56.3</td><td>16.4</td></tr><tr><td>Longitude</td><td>11.7</td><td>23.2</td></tr><tr><td>Depth</td><td>10.3</td><td>28.1</td></tr><tr><td rowspan="2">L. sarmentosa</td><td rowspan="2">0.9731</td><td>Longitude</td><td>38.7</td><td>65.6</td></tr><tr><td>Rugosity</td><td>5.6</td><td>10</td></tr></tbody></table>