Ecological role of benthic Vulnerable Marine Ecosystems (VMEs) indicator taxa on soft bottoms

Study area and data collection

The research was conducted in the northern sector of the Strait of Sicily, with the centroid of the study area located at 38°0.587’N, 11°19.329’E (Fig. 1). The area, part of the MedWind project framework, spans a total surface of approximately 1651 km². Acoustic mapping was performed between August 12th and September 6th, 2021, using two multibeam echosounders: a Kongsberg EM2040 for depths ranging from 150 to 300 m, and a Kongsberg EM712-MK2 for deeper zones reaching 1000 m. The resulting data were processed to generate a Digital Terrain Model (DTM) with a spatial resolution of 5 m.

ROV-based visual surveys were subsequently carried out between September 11th and November 17th 2021, over a 67-day period, from the MainportGeo research vessel. The remotely operated vehicle used was a Tomahawk Light Work Class ROV equipped with dual manipulators, four types of cameras (full HD, standard color and black-and-white, and a 6 K high-definition camera), laser scaling devices, a sampling box, beacon, DVL system, and a Seabird Microcat SBE 37. A total of 140 transects were completed across a depth range of 135 to 885 m, covering a total distance of 129.5 km (mean transect length: 929 m ± 257 SD).

During the video surveys, all observations were georeferenced and recorded using the Ocean Floor Observation Protocol (OFOP) software which logged time, date, ROV and vessel positions, depth, substrate characteristics, and encountered species. In parallel with the ROV deployments, oceanographic data (temperature and salinity) were collected at 97 randomly selected stations using a Rosette sampler and a Sea-Bird Scientific SBE 911 Plus V2 CTD probe. Of the 140 transects collected during the ROV surveys, a subset of 74 transects were selected for analysis based on their dominant soft-bottom substrate, which averaged approximately 90% coverage (Fig. 1 and Table S1).

Explanatory variables

The explanatory variables included in the models to assess the influence on the distribution patterns of associated fauna were classified into different groups: VME indicator taxa, morphological (slope, depth, rugosity, aspect), environmental (temperature, salinity and Chl-a), and anthropogenic (Bottom Trawling Fishing Effort – BTFE) (Figure S1).

The OFOP analysis shows that the selected transects were dominated by three VME indicator taxa (I. elongata, Pennatuloidea and L. phalangium). These taxa were counted (N), and their abundance standardized to one linear kilometer, calculating their density (N * Km⁻¹).

Morphological variables were derived from a high-resolution Multi-beam depth data (5 m) obtained through multibeam echosounder surveys. Slope, rugosity, and aspect were computed using the “terrain” function from the R package raster⁴³. Slope reflects the steepness of the seabed, ranging from 0° (flat) to 90° (vertical), and is known to influence benthic habitat distribution by enhancing local current flows and affecting fishing gear accessibility^44,45. Rugosity was calculated as the elevation difference between neighboring cells and provides a measure of seafloor complexity, with higher values indicating more rugged or rocky terrains. Aspect describes the orientation of the slope and is relevant for assessing exposure to prevailing currents⁴⁶.

Environmental variables were modeled using spatial interpolation techniques. Bottom temperature and salinity data were collected from 97 CTD stations during the survey period and interpolated via co-kriging, with bathymetry included as a co-variable to improve prediction accuracy¹⁶. Chlorophyll-a concentration data were retrieved from the Copernicus Marine Service and represent the annual mean for 2021, spatially matched to the study area.

BTFE was estimated from Automatic Identification System (AIS) data for 2021, processed following the methodology described by Russo et al. (2016)⁴⁷. The spatial distribution of bottom trawling activity was reconstructed by identifying fishing events through vessel speed and bathymetric filters, and cumulative Bottom Trawling Fishing Effort (BTFE) was calculated as the total number of fishing hours per 1 km² grid cell.

Response variables

The associated fauna observed during the video surveys, fish and crustacean species, was treated as response variables in the statistical models. For both fish and crustaceans, the sighting density distribution of each species was calculated (N * Km⁻¹). Next, the quartiles of these two distributions were calculated, and species falling below the first quartile (both fish and crustaceans) were removed from the analysis to minimize noise and improve the robustness of the results. However, for crustaceans, the limited number of species observed posed challenges in terms of data interpretation. To address this issue and reduce the differences in density between species (e.g., high abundance versus rare observed species), the final dataset for crustaceans was converted into a presence/absence format. This approach allowed equal weighting of all species in the analysis, ensuring a more balanced representation of the crustacean taxa.

Data analysis

To assess the relationship between associated fauna (fish and crustaceans), and explanatory variables, a Redundancy Analysis (RDA) was performed. This multivariate ordination method assumes a linear relationship between species responses and the ordination axes⁴⁸. In this study, each of the 74 transects was considered as a single analytical unit. For each transect, the mean values of environmental (temperature, salinity, chlorophyll-a), morphological (depth, slope, rugosity, aspect), and anthropogenic (BTFE) variables were calculated. Before proceeding with RDA, to reduce multicollinearity among explanatory variables and improve model performance, a Variance Inflation Factor (VIF) analysis was conducted. Variables exhibiting VIF values greater than 3 were considered highly correlated and were excluded from further analysis⁴⁹.

All analysis were conducted using the “vegan” package in R⁵⁰. To evaluate the contribution of each explanatory variable in shaping the associated faunal assemblage, we applied the envfit() function⁵⁰ to all predictor variables (VME, morphological, environmental, and anthropogenic). This approach projects each variable onto the RDA ordination and tests its statistical significance via permutation. Only variables with p-values < 0.05 were considered significantly associated with community structure and were retained for interpretation and graphical representation. The percentage of variance explained by each RDA axis was calculated using constrained eigen values⁵⁰. This method ensures that only the model-constrained portion of variance is reported, excluding unconstrained residual variation.

To determine which species of associated fauna were significantly associated with the RDA axes, we calculated Pearson correlations between each species and the RDA axes (RDA1 and RDA2), and then we used two-sided tests to check if the correlation between species and RDA axes were significantly different from zero, either positive or negative.

RDA biplots were generated using ggplot2⁵¹ and ggrepel⁵² R packages. Species scores, site scores, and centroids of explanatory variables were extracted with the scores() function. The RDA axes were annotated with the percentage of variance they explained to facilitate interpretation of the ordination.