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 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 (González-Wangüemert et al. 2018). 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 (González-Wangüemert et al. 2014, 2018). 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 macroalgal-dominated 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.

Methods

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).

Fig. 1
figure 1

Study area at Arrábida’s coast, in the southwest of Portugal (38°29′19″N 8°57′03″W), with sampling stations in red and the respective sampling station code

Full size image

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 (Félix et al. 2021). 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 (Félix et al. 2021).

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).

On each survey, for each station, water column parameters were recorded with a multiparametric sonde (YSI-EXO2, Xylem Inc., USA), namely temperature (°C), pH, salinity, dissolved oxygen (DO) (mg L−1), total dissolved solids (TDS) (mg L−1), turbidity (NFU), chlorophyll a (Chla) (μg L−1), depth (m) and longitude (°). Current intensity (m s−1) was measured using a Doppler Current Sensor (model 4100, model 4100 from Aanderaa, Norway). Sediment samples were collected to determine the total organic matter (TOM), obtained by loss of weight on ignition (480 °C) (Byers et al. 1978). Sediment granulometric fractions were determined using 63-μm to 2.0-mm sieves to separate the silt, sand and gravel fractions. Each fraction was dried and weighed. Sediments were ranked according to their percentages (Blott and Pye 2001) and classified using Shepard diagrams (Shepard 1954).

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 m2, 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 (Félix et al. 2021). 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).

Holothuria arguinensis showed an irregular mean density distribution through time, but with no statistically significant differences. There was, however, an increase in density in January, March, and August (Fig. 2a). The average observed density for the distribution range was 0.69 ind. 100 m−2, and the maximum registered density was 10 ind. 100 m−2 at SE10 (Fig. 2b). Mean density showed statistically significant spatial differences (KW: p = 4.352e−7) with both SC2 and SE10 having the highest densities (1.18 and 1.59 ind. m−2, respectively) (Fig. 2b). The SE10 station differed significantly from all other stations (Dunn: p < 0.05) except SC2 (Dunn: p = 1), with the most different pair being SE10 and SC3 (Dunn: p = 1.6045e−05). There was also an irregular spatial mean size distribution (ANOVA: p = 0.0057) with SC3 and SE10 having the biggest individuals (Fig. 3c). SC6 and SE10 were the most different stations when comparing sizes (Tukey: p = 0.00190) followed by SC3–SC6 (Tukey: p = 0.04297). However, there were no statistically significant size differences over time (ANOVA: p = 0.437) (Fig. 3).

Fig. 2
figure 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

Full size image
Fig. 3
figure 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

Full size image

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 R2 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.

Table 1 Final GLM results with a Tweedie distribution (R2 = 33%) explaining density distribution patterns in the species Holothuria arguinensis, in the Arrábida coastal area, Setúbal, Portugal
Full size table

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.

Fig. 4
figure 4

Surface plots, representing the behaviour of the interaction between a substrate and depth and b current and depth in the final model, explaining the density distribution patterns of Holothuria arguinensis in Arrábida. The substrate is a binary variable, with 0 representing the rocky substrate and 1 the sandy substrate

Full size image

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).

Table 2 Full average model results with ∆AIC < 2 (N = 11), explaining the size distribution patterns of Holothuria arguinensis in the Arrábida distribution range, Setúbal
Full size table
Table 3 Values of importance of variables in the final full averaged model used to explain the size distribution of Holothuria arguinensis, considering the 11 main models with a ∆AIC < 2
Full size table

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 (Félix et al. 2021) 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 (Félix et al. 2021).

There was a clear preference for rock transects but when considering the actual surfaces, regardless of transect, H. 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 (Félix et al. 2021), 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 (Félix et al. 2021). Holothuria arguinensis does not require shelter in the same way as H. mammata because it grows faster in weight with volume, potentially avoiding 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 (Félix et al. 2021), 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 (Félix et al. 2021).