Spatiotemporal Distribution Patterns of the Sea Cucumber Holothuria arguinensis on a Rocky-Reef Coast (Northeast Atlantic)

The demand for sea cucumbers from the NE-Atlantic is increasing as a consequence of the depletion of stocks that historically supplied Asian markets. Some species from the NE-Atlantic are already showing overfishing signs due to either inefficient policies or their total absence, and they do not offer any margins of sustainability to sea cucumber fisheries. The current lack of knowledge of biological and ecological traits has led to ineffective protection of those species. Holothuria arguinensis is a good example of a temperate sea cucumber species suffering from the effects of unregulated and illegal harvest. This study aims to determine the temporal and spatial distribution patterns of H. arguinensis, assessing density and size-class structure in a NE-Atlantic population, in SW Portugal, as a function of environmental conditions. Sampling was done every 1.5 months at 6 stations defined according to habitat heterogeneity, following the coastline from the mouth of the Sado estuary to marine coastal areas. At each station, there were two habitats (rocky substrate and sandy areas) sampled through SCUBA diving along 30-m random transects, each with three replicates. Individuals were counted and measured (total length) and individual behavioural parameters, such as sheltering and settlement surface, were recorded. Water column parameters, sediment composition, and substrate cover were also recorded. Generalised Linear Models were fit to the data, modelling density and size distribution of the species with respect to environmental conditions. The results indicate that the distribution of H. arguinensis was mainly determined by its ability to deal with physical environmental stressors (current), rather than the chemical properties of the water column or sediment composition. The hydrodynamic features of the region exert the greatest impact on this species’ distribution and help shape habitat usage by different size class groups.


Introduction
Sea cucumbers are a diverse group of echinoderms with a key ecological role in marine ecosystems. Deposit-feeding sea cucumbers promote bioturbation through their feeding Communicated by Linda Deegan * F. Azevedo e Silva fhsilva@fc.ul.pt activities, recycling and distributing nutrients, contributing to the health of benthic habitats and, consequently, its biota (Uthicke 1999;MacTavish et al. 2012;Purcell et al. 2016). However, as a consequence of commercial fishing and despite their positive effects on surrounding marine communities, several sea cucumber species have decreased in numbers, potentially creating cascading negative effects in the ecosystems. This reduction in sea cucumber populations is a consequence of overfishing driven by their high value in Asian markets and increasing demand (Máñez and Ferse 2010;Purcell et al. 2013;González-Wangüemert et al. 2018). The overexploitation of several sea cucumber species, particularly in the Indo-Pacific, has led to an expansion of fisheries to new areas, targeting new species (Purcell et al. 2013;Eriksson and Clarke 2015;González-Wangüemert et al. 2018). Holothuria (Roweothuria) arguinensis is a good example of this, having achieved high economic value and being intensively and illegally harvested in the south of Portugal and Spain ). In some regions, a decrease in density, abundance and genetic diversity, as well as loss of the largest size classes, increase in prevalence of some diseases and local stock depletion have already been linked to the fishing pressure . Because of recent economic interest in this species, some studies have addressed the population status, considering both its reproductive cycle (Domínguez-Godino et al. 2015;Marquet et al. 2017) and its ecological traits (González-Wangüemert et al. 2013;González-Wangüemert and Borrero-Pérez 2014;Navarro et al. 2014;Siegenthaler et al. 2015;Domínguez-Godino and González-Wangüemert 2020). However, these studies all focused on intertidal populations at lower latitudes. Given the species plasticity (Hadfield and Strathmann 1996), environmental variations can induce physiological adaptations and, consequently, changes in life-history traits. Holothuria arguinensis is a temperate sea cucumber distributed along the Northeastern Atlantic from Berlengas Island (West Portugal) (Rodrigues 2012) to Morocco and Mauritania, including the Canary Islands (Tuya et al. 2006;González-Wangüemert and Borrero-Pérez 2014) and western Mediterranean Sea (González-Wangüemert and Borrero-Pérez 2014). In those locations, this species occurs more frequently in the intertidal zone, on macroalgaldominated beds and in seagrass meadows where individuals often find cover (Navarro 2012;Domínguez-Godino and González-Wangüemert 2020). Holothuria arguinensis also reveals density-and size-dependent behaviour where bigger individuals appear linked to the previously described habitats and with a preference for the lower intertidal zone (2 m deep at low tide) and smaller individuals in the upper intertidal (1 m deep at low tide) (Domínguez-Godino and González-Wangüemert 2020). Seasonal behaviour was also described for intertidal areas, where higher environmental variability leads individuals to migrate to deeper water in the summer, to avoid excessive heat and UV radiation exposure (Domínguez-Godino and González-Wangüemert 2020). However, there is a gap in knowledge on how this species uses non-intertidal habitats and different habitats, such as rocky reefs or sandy areas, and how they might induce different behaviours.
The aim of this study was to determine the temporal and spatial patterns of H. arguinensis density and size distribution in a NE-Atlantic coastal subtidal area, in SW Portugal, as a function of environmental characteristics. Size and density distributions were modelled to describe its habitat preferences.

Study Area
The study took place in Portugal, in the district of Setúbal, in the coastal area adjacent to the Sado estuary, more specifically in the Professor Luiz Saldanha Marine Park, a Marine Protected Area (MPA) (38°26′50.4′′N; 9°01′58.7′′W), and near the estuary mouth (Fig. 1). Sampling stations were selected according to habitat heterogeneity, determined by the biophysical and geomorphological characteristics of the site. They were also selected according to the species distribution area, as defined by Félix et al. (2021), which assessed the distribution area of commercial sea cucumbers in different habitats and different depths. A preliminary exploratory survey, in addition to the information gathered from Félix et al. (2021), allowed identification of locations (stations) where H. arguinensis occurred which constitutes the spatial extent of the present study. For this study, the distribution range was restricted to rocky reefs and adjacent areas. The coastal area south of the estuary, including phanerogamic grasslands, is exclusively composed of soft substrate and has no sea cucumbers. The stations where H. arguinensis occurred were SC2, SC3, SC4, SC5 and SC6 (marine stations), with only one estuarine station (SE10) (Fig. 1).
The Sado estuary is more heterogeneous in environmental features than the adjacent coast. Dissimilarities from the adjacent marine environment were mainly related to the granulometric characteristics of the sediment, sediment organic content, turbidity, dissolved oxygen and temperature . The only location with sea cucumbers within the estuarine area was SE10, a rocky outcrop in the estuary mouth. This rocky outcrop resembles the Arrábida coastal rocky reefs, unlike the remaining estuarine areas .
The marine area of the study is inside the MPA. The Arrábida MPA has particular characteristics contributing to high biodiversity. First, it is near the northern limit of a major upwelling system near the Sado estuary, which is the source of increased phytoplankton biomass during the spring and summer (Santos et al. 2022). It is mostly composed of rocky reefs with a gradual transition to sand bottoms. Finally, the area mostly faces south, so it is sheltered by the Arrábida mountain range from prevailing winds from the north and from the northwest swell, the most frequent swell orientation on the Portuguese west coast (Wooster et al. 1976;Costa et al. 2013). Considering differences between marine habitats, six sampling stations were established based on hydrodynamic characteristics (current intensity), availability of shelter (e.g. crevices in rocks), and distance to the estuary.

Sampling
Stations were surveyed every 1.5 months from Jan 2018 to Mar 2019. Counts per unit area were recorded using the same sampling method at all stations. SCUBA diving surveys were done between 09:00 and 12:00 (visual census in randomly placed 30 × 30-m transects, parallel to the coast), with three replicates in each sandy and rocky substrate. Each individual was carefully measured in situ (total length) with a flexible measuring tape, avoiding handling-related stress that can induce muscle contraction and change a sea cucumber's length (Azevedo e Silva et al. 2021). The surface where each individual occurred was noted, i.e. a sea cucumber on a rocky substrate transect could be on a sandy patch or on rock and vice versa. The presence of shelter and cryptic behaviour of each sea cucumber were also recorded (e.g. macroalgal cover or crevices and chambers in rocks).

Demographic and Density Distribution Data Analysis
Data analysis was in two phases. First is an exploratory identification of space, time and substrate differences in the density and length distribution of the population. This allowed for the early identification of key patterns in the species distribution in the area. Then, the distribution patterns of the species were evaluated in terms of density and total length (a proxy for demography) as a function of the environmental conditions in the distribution range, via Generalised Linear Models (GLM).
In the first phase, the existence of spatial, temporal and substrate differences in the density distribution (non-normal distribution, with positively skewed data) was evaluated with a Kruskal-Wallis test. A posteriori paired multiple comparisons were performed using the Hochberg method (with an adjusted p value) (Hochberg 1988). Temporal and spatial differences in the size (normally distributed, Shapiro-Wilk: W = 0.97872, p value = 0.1422) according to substrate or cryptic behaviour were evaluated by an ANOVA, followed by a posteriori Tukey analysis if required.
In the second phase, to determine habitat preferences of the species as a function of environmental conditions, Generalised Linear Models were fitted to the data. To avoid misleading inferences due to similar effects of correlated variables in a smaller sample size data set, a selection was made based on Variance Inflation Factors (VIF), with a threshold set at 8 and then checked with correlation plots. Fine and medium fractions of sand and gravel were excluded since coarse sand and silt best explained granulometry of the sediment. Both salinity and total dissolved solids (TDS) were highly correlated and were excluded, because first salinity in coastal areas is less variable than in estuarine environments and salinity at the SE10 site, the only estuarine station, is influenced by the tide. Then, sea cucumbers have low mobility, and since sampling could not be standardised according to tide, salinity and TDS variations, which are tide dependent, were poor predictors for distribution patterns. Other variables such as longitude, turbidity, and dissolved oxygen were chosen to represent distance to the estuary. Longitude represents the most reliable measure of the potential influence of the estuary (based on distance), since it does not change, unlike the physical-chemical parameters of the water column which are influenced by season and tide. Shelter represents the presence of crevices or seams, being a proxy for a behavioural response to stressors such as current or predation. The local presence of shelters was strongly correlated with the variable current, causing convergence issues in the models, as the sites without shelter were also those with the greatest exposure to current influence. Thus, the variable shelter was removed, since current was a better predictor.
For these analyses, GLMs were implemented for two response variables: density and mean total length (both per replicate). The residuals of mean lengths were not normally distributed, and for this variable, the analyses were fitted with Gamma distribution and a logarithmic link function. The density distributions showed a zero-inflated distribution, with different CPUE records depicting a point mass at zero, a positive tail and non-negative values. This led to the choice of a Tweedie distribution-that can deal with excessive zeroes-for the response (Tweedie 1984;Shono 2008). For both variables, after fitting the full models, model selection was made through an evaluation of all possible combinations of variables via the dredge() function from the package MuMIn (Bartoń 2020) and based on the ∆AIC (AIC (Akaike Information Criterion)). Models with a ∆AIC < 2 were considered for further inference, as having similar explanatory power. If more than one model presented ∆AIC < 2, then averaged inferences via model.avg() provided the model-averaged coefficients (full and conditional averaging). Since model averaging makes interpretation using p values to test a particular variable more difficult, they were additionally subjected to interpretation based on the weight of the predictors in the average model (Grueber et al. 2011). The full model coefficients set terms to zero when they are not included in a given model whilst averaging, whereas the conditional coefficients ignore the predictors whenever they are not included in a model and only consider them in the models where they are represented. Thus, the full model coefficients are more conservative (Burnham and Anderson 2002). The interactions between significant variables in the models were also analysed, determining their type, considering the signs of the coefficients, and the behaviour of one predictor as a function of the other, using surface plots (Feld et al. 2016). All statistical analyses were performed using the R software (R Core Team 2020).

Results
A total of 174 individuals were recorded, in a total sampled area of 29,160 m 2 , with lengths of 80−490 mm. Only three individuals were found in the study area in habitats with macroalgal cover. Of all the individuals, only 3% were found sheltered at the time of sampling (daytime). Most individuals were found in rock transects (90%) as opposed to sand transects (10%) (Kruskal-Wallis: p = 1.179e −11 ). However, if we consider the individual settlement surface, which is independent of transect, there was no clear preference for sandy or rocky substrate at the time of sampling (50%/50%). It is, however, important to acknowledge that sand transects had few rock outcrops, and in rock transects, sandy patches were common . When considering size distribution between settlement surfaces, individuals found on rock were not shorter than individuals on sandy patches (ANOVA: p = 0.285). However, considering habitats, sea cucumbers on rocky transects were shorter (303.1 mm ± 6.7) than individuals on sandy transects (363.9 mm ± 10.7) (ANOVA: p = 0.000339).

Density Distribution Models
To explain the density distribution of H. arguinensis and after running all possible combinations of the predictive variables, the GLM simplification resulted in a single model with an ∆AIC > 2, an adjusted R 2 of 33% and an explained deviance of 39.6% (Table 1). Significant variables in the fitted model were (1) depth, (2) substrate, (3) current, and (4) % silt. Density, lower in sand transects, decreased with increasing depth, decreased with current and decreased with increasing percentages of silt in the sediment.
There were also two significant interactions in the final model, both with depth: depth and substrate; depth and current. Both interactions were antagonistic, meaning an increase in depth cancelled the negative effects of both current and sand substrate on density. The first interaction revealed the decrease in density in sand transects was less pronounced with an increase in depth (Fig. 4a). The second interaction revealed that an increase in depth cancelled the negative effect that current had on density (Fig. 4b), i.e. at higher depths density did not decrease with increasing current, unlike in shallower areas.

Size Distribution Models
Based on the ∆AIC < 2 criterion, 11 models were selected to explain the size distribution of H. arguinensis. Because full and conditional averages of these models lead to similar results, only the results of the full average of the models are presented since it is more conservative ( Table 2). The most significant variables in model fitting were coarse sand (%), with larger individuals on coarser sand; substrate (sand), with larger individuals on sand transects; and longitude (distance to estuary), where larger individuals were closer to the estuary. When considering the 11 models, the sum of weights of the variables revealed only coarse sand (%), substrate (sand) and longitude in all models (Table 3). Between the remaining and less contributory variables, the organic matter of the sediment was the most relevant (values up to a maximum of 8.5% TOM).

Discussion
There are three sea cucumber species in the studied area, Holothuria arguinensis, Holothuria mammata and Holothuria forskali (Azevedo e Silva et al. 2021; Félix et al. 2021). In the present study, the densities of H. arguinensis were lower than in studies with the other two sympatric species in the same area. In fact, the densities of H. mammata ) and H. forskali (unpubl. data) were > 10 times than what was reported in this study, for H. arguinensis. The maximum estimated density of H. arguinensis (10 ind. 100 m −2 ) was at the rock outcrop at the estuary mouth. This density is higher than what was found in other studies in the south of Portugal (maximum average of 6.8 ind. 100 m −2 at Ria Formosa) (Siegenthaler et al. 2015;Domínguez-Godino and González-Wangüemert 2020) and slightly less than in the Canary Islands (maximum of 12 ind. 100 m −2 ) (Navarro et al. 2013).
With a maximum total length of 490 mm, the present study area had the largest individuals of this species. Other studies in the Mediterranean and Northeastern Atlantic recorded individuals 70−410 mm, with the previous highest record in the south of Portugal (Siegenthaler et al. 2015;Marquet et al. 2017;Olaya-Restrepo et al. 2018;Domínguez-Godino and González-Wangüemert 2020). The population from Ria Formosa has been subjected to illegal harvesting (Olaya-Restrepo et al. 2018) which typically results in the loss of the bigger size classes. Considering that in the present study, H. arguinensis is at its northern distribution limit, the evident differences in sizes and the densities found are a strong indicator that, in the study area, this is a healthy, unexploited population. Unlike in the areas from other studies, access to the Arrábida coast is difficult. The area also benefits from habitat complexity favoured by upwelling events in the summer and nutrient inputs from the estuary to the marine environment throughout the year (Wooster et al. 1976;Cabeçadas et al. 1999;Costa et al. 2013;Santos et al. 2022). The difficult access to the coast and the upwelling events prevent illegal harvesting and generate high productivity, respectively, which contribute to the good condition of the population. Also, most previous studies describe occurrences in intertidal zones, whilst in this study area this species occurs exclusively in subtidal habitats. Nutrient-rich subtidal environments provide more stable conditions for sea cucumber populations and this could be the main reason for higher abundances and bigger sizes with lower temporal variation, as described for H. mammata ).
There was a clear preference for rock transects but when considering the actual surfaces, regardless of transect, H. Fig. 2 Densities of Holothuria arguinensis for its distribution range in the Arrábida region, during the sampling period (a) and by sampling site (b). In (a), the red dots connected by lines represent the evolution of the mean, and in (b) * represent the site mean 1 3 Fig. 3 Total lengths registered for Holothuria arguinensis, for its distribution range in the Arrábida region, for the total area (a), during the sampling period (b) and by sampling site (c). In (b), the red dots connected by lines represent the evolution of the mean, and in (c) * represent the site mean 1 3 arguinensis was found equally on sandy patches or rocks. Sandy patches within rocky transects always had shelter close by, as rocky transects represented a more complex substrate, often with crevices and chambers. Sand transects were open and exposed with little or no shelter nearby. The differences in sea cucumber size found between habitats (rock or sand transects) suggest a life strategy that is dependent on size. Migration between habitats, from juvenile to adult, has been described for other species such as Holothuria scabra. Juveniles of this species stay on leaves of seagrass and then migrate to sandy areas as adults (Mercier et al. 2000). Recently, for H. arguinensis, a size-dependent strategy has also been described, where smaller individuals preferred a more complex habitat with different microhabitats across small areas and the bigger individuals remained in more exposed areas (Domínguez-Godino and González-Wangüemert 2020). Considering that in the study area, H. arguinensis feeds almost exclusively on sand (Azevedo e Silva et al. 2021), this could be an evolutionary response to predation and environmental stressors where larger sea cucumbers move opportunistically into open spaces to feed. Smaller sea cucumbers stay in more complex habitats, near shelter, as the risk of predation and negative impact from environmental stressors decreases with size (Shiell and Knott 2008;Purcell 2010;Ceccarelli et al. 2018). This is similar to the behaviour of H. mammata in the same area ), but in different proportions, as H. mammata is more often found sheltered than H. arguinensis. Holothuria mammata habitat distribution preferences include environmental stability and presence of shelter where smaller individuals occupy more complex habitats and bigger individuals venture into open, more exposed areas to feed . Holothuria arguinensis does not require shelter in the same way as H. mammata because it grows faster in weight with volume, potentially avoiding  The substrate is a binary variable, with 0 representing the rocky substrate and 1 the sandy substrate predators and high currents more efficiently (Purcell and Simutoga 2008;Azevedo e Silva et al. 2021). Such behaviour (no shelter) has also been reported for H. arguinensis in the Canary Islands (Navarro et al. 2014).
The GLM results revealed that abiotic and biotic parameters of the water column were not the main factors driving H. arguinensis abundance and size distribution, as reported by Domínguez-Godino and González-Wangüemert (2020) for a population in southern Portugal and by Félix et al. (2021) for H. mammata, in Arrábida. Depth and current were the most influential factors affecting the abundance of H. arguinensis. Like one other species in the area (H. mammata; Félix et al. 2021), H. arguinensis is affected by higher currents, preferring a more stable environment. This is noticeable in the antagonistic relationship between the variables depth and current, where the negative effect of the physical disturbance by the current was cancelled by an increase in depth, suggesting that, in areas where currents are stronger, individuals may seek refuge in deeper waters. The highest densities found at the rock outcrop at the mouth of the estuary were a good reflection of this behaviour, as this location is a place of swift currents (influence of the tide on the estuary bar) and also one of the deepest sites in this study. A high phytoplankton biomass, associated with the input of nutrients from the Sado estuary (Cabeçadas et al. 1999;Santos et al. 2022), may also favour feeding conditions and promote greater aggregation of individuals of this species at greater depths.
The habitat usage strategy of this species shows some differences from the other species in the study area , notably the density distribution of H. arguinensis is not related to size. Higher densities do not represent younger groups, but they suggest an opportunistic strategy that could be related to better feeding conditions. Results indicate that bigger individuals are found closer to the estuary, in sandy habitats, and in locations where the sediment is coarser, indicating that different sizes have different habitat use strategies. Even though individuals of this species occur mainly on the rocky reef, the larger individuals were more frequently found either in stretches of sand within the reefs or even in sandy habitats, more exposed. The SC3 station, the least favourable environment when considering the most influential variables on density distribution, had the biggest individuals. Larger individuals are less likely to be impacted by tidal currents and less likely to be preyed upon, which gives them the possibility of occupying areas with higher currents.
Despite the preference of this species for sediments with a higher percentage of coarse sand, the percentage of organic matter present in the sediment may also indicate a preference for sediments with higher organic content, even though the percentage of organic matter in the sediment is nowhere near as impactful as the hydrological characteristics of a site. This tendency has also been described in habitat selection by this species at Ria Formosa (Domínguez-Godino and González-Wangüemert 2020) and in the Canary Islands (Navarro et al. 2014) and even in other detritivores (Slater et al. 2011).
In general terms, the size and density distribution of H. arguinensis is clearly shaped by its ability to deal with physical environmental stressors, rather than the biotic and chemical properties of the water and substrate. This allows it to occur in sympatry with other species, avoiding competition for space and resources. It also allows for size-dependent strategies, culminating in different habitat selection by individuals of different sizes. This could increase the fitness of the species by expanding available habitat for resource exploitation. When compared to previous research, it is also important to note that, despite being conducted in a constrained area on a North Atlantic reef, there is consistency in the factors that influenced the distribution of various detritivorous sea cucumbers (e.g. Sloan and von Bodungen 1980;Sonnenholzner 2003;Mendes et al. 2006;Morgan 2011;Domínguez-Godino and González-Wangüemert 2020;Félix et al. 2021). These factors, along with details on the species' growth and reproductive cycle, will provide guidance for the future management of the species, as was noted for H. mammata ).
Funding Open access funding provided by FCT|FCCN (b-on). This work was funded by the Operational Programme Mar2020 MAR-02.01.01-FEAMP-0052 "Newcumber-Avanços para o cultivo sustentável de pepinos do mar". This work was also supported by national funds through FCT-Fundação para a Ciência e a Tecnologia, I.P., within the strategic projects MARE (UIDB/04292/2020) Associate Laboratory

Data Availability
The data used in this study are available by request from the corresponding author (F. Azevedo e Silva).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.