Modelling the Distribution of a Commercial NE-Atlantic Sea Cucumber, Holothuria mammata: Demographic and Abundance Spatio-Temporal Patterns

Introduction

Sea cucumbers are a group of echinoderms with important ecological functions. Deposit-feeders compose a large group of these holothurians, which play a critical role in habitat structuring through bioturbation, recycling and redistributing nutrients (e.g., MacTavish et al., 2012; Purcell et al., 2016). Their feeding behaviour, ingesting sediments and digesting its organic content, help to maintain or improve the physico-chemical processes of benthic habitats, by decomposing sediment organic matter and controlling populations of bacteria, fungi, and phytobenthos (e.g., Uthicke, 1999; Michio et al., 2003; Namukose et al., 2016). Sea cucumbers are a highly priced commodity in China and have been harvested for more than three centuries (Schwerdtner Máñez and Ferse, 2010). As a consequence of high fishing pressure and intrinsic biological and ecological characteristics, such as low mobility, density dependent breeding, slow growth of many species, or shallow water occurrence for others (Hamel and Mercier, 1996a; Kinch et al., 2008; Mercier and Hamel, 2009; Purcell et al., 2013), this group is vulnerable to overfishing. The majority of exploited sea cucumber species are, in fact, deposit-feeders (most of the formerly known aspidochirotids) (Purcell et al., 2016), and several of these are currently overexploited or depleted (Purcell et al., 2013). This can create a set of inevitable ecological consequences that can result in a curtailment in primary productivity, decrease in benthic diversity and biomass, decrease of buffering capacity of surrounding seawater acidity, preclude symbiotic associations, impair trophic relationships, etc. (Moriarty et al., 1985; Uthicke, 1999; Michio et al., 2003; Solan, 2004; So et al., 2010; Wolkenhauer et al., 2010; Schneider et al., 2011, 2013; Purcell et al., 2016). Considering the ecosystem services provided by these species, the overexploitation of multiple populations has the potential to create severe cascading effects in the ecosystems.

Setting limits to fisheries is now urgent but, for that to be done efficiently, improved knowledge on the ecological traits of exploited species is of the utmost importance. Distribution patterns and their drivers are a fundamental set of information for stock management, e.g., establishing geographical limits to harvest and avoid excessive capture of breeders or specific size-classes, particularly due to a common patchy distribution of sea cucumbers, where high density areas can be followed by adjacent zero or near-zero abundances (e.g., Džeroski and Drumm, 2003; Mendes et al., 2006; Eckert, 2007; Shears and Babcock, 2007; Dissanayake and Stefansson, 2010; Domínguez-Godino and González-Wangüemert, 2020). The study of these traits takes on a new importance considering that this is not only species-specific but varies according to habitat, i.e., the distribution patterns of a species may depend on the habitat it occurs in Sloan and von Bodungen (1980), Navarro et al. (2013), Domínguez-Godino and González-Wangüemert (2020). Depending on habitat type, densities of deposit-feeders may be higher in shallow-waters, but may vary in depth, for reasons like reduced exposure from UV radiation (Domínguez-Godino and González-Wangüemert, 2020). A vegetated cover can be the main driver for aggregation, as it provides shelter (Dissanayake and Stefansson, 2012; Navarro et al., 2014; Siegenthaler et al., 2015, 2017) or the opposite (Tuya et al., 2006), but in some cases the preference can change according to life-cycle stage (Mercier et al., 2000; Hamel et al., 2001; Eriksson et al., 2012; Aydin, 2019a). Rocky reefs can also be preferred when compared to adjacent areas without rocky substrate, offering more complex habitats, regarding food sources or shelter (Zhou and Shirley, 1996; Džeroski and Drumm, 2003; Mendes et al., 2006). Generally, density is not related to organic content of the sediment, although it can be related to its granulometry and, consequently, to diet quality through benthic biodiversity (Sloan and von Bodungen, 1980; Mercier et al., 2000; Dissanayake and Stefansson, 2012; Navarro et al., 2014; Domínguez-Godino and González-Wangüemert, 2020; Viyakarn et al., 2020). Muddy sediments, with more organic content, have, typically, less microbenthic diversity than sandy sea beds (Underwood and Barnett, 2006). Often, the main habitat preferences stem on lower hydrodynamic characteristics of the area and the presence of shelter (Sloan and von Bodungen, 1980; Mendes et al., 2006; Al-Rashdi et al., 2007; Morgan, 2011).

However, and in spite of the ecological relevance of these species and risk of overexploitation, these traits are fairly described for the most relevant species in the trade market (e.g., Mercier et al., 2000; Shiell, 2007; Morgan, 2011; Dissanayake and Stefansson, 2012; Purcell et al., 2012, 2016; Kashio et al., 2016) but for others are sparse and scarce. Emergent target species (from the Mediterranean and NE-Atlantic) are less studied, or not at all, or have a strong regional or habitat focus (e.g., Simunovic et al., 2000; Navarro et al., 2014; Siegenthaler et al., 2015, 2017; Marquet et al., 2017; Aydin, 2019b; Boncagni et al., 2019; Domínguez-Godino and González-Wangüemert, 2020). The introduction of these species in the international trade markets is the result of a stock depletion of traditional commercial species, mainly from the Indo-Pacific, increasing demand and consequent fisheries expansion to new areas (Purcell et al., 2013; Conand, 2017; González-Wangüemert et al., 2018; Dereli and Aydın, 2021). Holothuria mammata Grube, 1840, is currently one of the most relevant commercial species in the NE-Atlantic and Mediterranean (González-Wangüemert et al., 2014, 2018; Dereli and Aydın, 2021). Its geographical distribution covers the NE-Atlantic where it has its northern limits at the Bay of Biscay and southern limits at the meridional Macaronesian islands including the insular regions of Madeira, Azores and Canary islands, the north-African coast, occurring also in the Mediterranean sea (Tortonese, 1965; Borrero-Pérez et al., 2009). Still, the knowledge on ecological traits for this species is scarce, limited and regional, and insufficient to establish stock management policies, although relevant efforts have been put into it (Navarro et al., 2013; González-Wangüemert et al., 2014, 2016, 2018; Aydin and Erkan, 2015; Mustapha and Hattour, 2017; Siegenthaler et al., 2017; Azevedo E Silva et al., 2018; Félix et al., 2018; Simões et al., 2018; Aydin, 2019a,b, 2020; Sousa et al., 2020; Venâncio et al., 2021; Azevedo e Silva et al., 2021). There are several ecological details missing, either to allow management and control of harvest activities or for future mitigation actions involving population recovery or replenishment. Hence, there is a need for further studies on reproductive parameters (across the species distribution range), growth, recruitment, population structure and dynamics, meta-population variability, and other life-history traits.

This study aims to determine the temporal and spatial patterns of distribution of H. mammata, considering its abundance and demographic structure in a NE-Atlantic area, SW Portugal. Individual size and density will be modelled as a function of the site’s environmental characteristics to describe habitat preferences.

Materials and Methodology

Study Area

The study was carried out at a NE-Atlantic coastal area in the southwest of Portugal, within a marine protected area (MPA), the Arrábida Marine Park (38°26′50.4″N; 9°01′58.7″W). The study area is dominated by a rocky reef that ends in a sandy sea floor in a gradual transition around 3–10 m in depth, depending on site, as depth varies greatly throughout this coastline. This subtidal area displays a variety of habitats in complex diversity of macro- and microhabitats, which supports a high diversity of algae, invertebrates and fish that benefit from the hydro- and geomorphological conditions of the area. The productivity of these coastal waters increases in the summer, as a consequence of upwelling events (Wooster et al., 1976; Costa et al., 2013). Sampling was extended offshore up to the 20-meter bathymetric mark and the Sado estuary to delimit the actual range of occurrence of H. mammata (Figure 1). This is a currently unexploited population. Fisheries are mostly banned in this MPA, where local marine authorities display a regular surveillance for illegal activities, and sea cucumbers, in particular, are not yet a fishing target.

Figure 1. Map of the study area, with sampling sites, grouped according to habitat type, classified as; estuary (yellow); coastal (red); and offshore (blue).

Sampling

Considering the lack of knowledge on the distribution of H. mammata extensive preliminary surveys were conducted considering habitat type, depth and area to identify the distribution of the species in the region (Figure 1 – all sites). These surveys, that allowed to determine which areas should be considered for the bulk study, continued throughout a 1-year period, at each of its four seasons (in January, May, August, and November), to account for temporal variability and guarantee the sites with no occurrences were true zeroes. The majority of estuarine sites go from muddy to sandy sediment, without cover, except for sites SE7 and SE8 that were set within subtidal patches of Zostera marina, Linnaeus. Sites SE9 and SE10 were also exceptions regarding the spatial distribution of sampling sites. They are at a close distance but have markedly different properties: SE10 is the only site within the estuarine area with a subtidal isolated rocky outcrop that resembles the coastal rocky reefs; and SE9 is an open area without shelter or biotic cover of the substrate. The marine area (offshore and southern coast) is dominated by a sandy substrate regardless of depth. Sites SC2 to SC6 differ from the remaining sites by the presence of a rocky reef (Supplementary Table 1).

Based on these surveys, each site with a recorded presence of H. mammata was included in the distribution range of the species. This distribution area included the sites SC2 to SC6 and SE10, which were then specifically set based on spatial heterogeneity considering hydromorphological characteristics (exposure to tidal currents), available substrates (sand, rock, density of algae cover in each of the previous), shelter frequency (available crevices in rock and in sand/rock transition areas) species densities and distance to the estuary. As this area of occurrence was the one used to assess the distribution patterns of the species, sampling effort was increased to 1.5-month intervals, between the months of January 2018 and March 2019 (1.5 sampling intervals were set to match the seasonal samplings that occur every 3 months). Thus, for sites SC2 to SC6 and SE10, 9 sampling campaigns were performed, while for all other sites, only four sampling campaigns were done, in total.

To ensure comparability across the study area sampling campaigns were always based on a catch per unit of effort (cpue) of number of individuals per unit of area, always carried out between 9 to 12 h (a.m.) and with three replicates per site. Sampling was primarily done by scuba-diving (visual census on random transects of 30 m × 3 m, parallell to the coastline), except for sites at depths higher than 12 m, with high currents or in upper estuarine areas, where visibility did not allow a visual census. In these sites, sampling was carried by beam-trawl with an opening of 2 m and a 2 mm mesh size at the cod end. Depending on location, tows of 5–10 min, also with three replicates per site, were performed. Typically, longer tows were possible offshore (higher depths and sandy bottoms) and shorter tows inside the estuary, with muddy bottoms and debris that promote a quick trawl clogging (Supplementary Table 1). Sites SC2 to SC6 and SE10 (species distribution area), unlike every other site, had two bottom types (rock and sand). For these, six random transects were done (three per substrate type) to represent the full heterogeneity of the site. The first transects were based on the deployment location and the following to the left and right of that location, parallell to the coastline. Transect areas were not repeated along the study period.

All H. mammata individuals detected in each transect were counted, and measured (nearest mm ±1), in situ, without manipulation to prevent muscle contraction and minimise size variability (using a flexible measuring tape). Since rocky substrates have sandy patches and vice-versa, the surface on which each individual occurred was also recorded (e.g., presence of an individual in sand substrate/patch of a rock transect), as was the cryptic condition of each individual (within macroalgal cover or sheltered in crevices and seams in rocks).

For each site (coastal, estuarine, and offshore) and each survey several environmental parameters of interest were measured, using a calibrated multiparametric sonde YSI-EXO2: water temperature, pH, salinity, dissolved oxygen (ODO), dissolved solids, turbidity, chlorophyll a, and depth. Site current was measured with a Doppler Current Sensor 4100. The type of rock surface for each site (plain or presenting crevices and seams) and the type of biological cover of the substrate (seagrass, macroalgae, or none) were also recorded. Grab samples were collected at each site for determination of sediment grain size and total organic matter (TOM) of sediments. The later was obtained by loss on ignition (480°C) and the former determined using 63 μm to 2 mm sieves to separate, respectively, the silt, sand and gravel fractions. Each fraction was then dried and weighed, and sediments assorted according to their percentages (Blott and Pye, 2001) and then classified with Shepard diagrams (Shepard, 1954).

Data Analysis

The data analysis was divided in two stages: (1) the exploratory assessment of occurrence areas for H. mammata and the identification of space and time differences in the distribution of the population; and (2) the use of empirical models to evaluate which environmental factors influence population density and size distribution of H. mammata. All statistical analyses were carried out in R software (R Core Team, 2020) and results considered significant at a p-value < 0.05.

To assess the environmental characteristics of the study area and to help explain the distribution range of H. mammata, a Principal Components Analysis (PCA) was used to group sites by their environmental parameters, considering three main habitats, classified as: estuary, coast, and offshore. The existence of spatial and temporal differences in the distribution of H. mammata, based on density and mean size, was assessed with a Kruskal–Wallis rank test. The selection of non-parametric techniques was justified by a non-normal distribution of density and a positively skewed distribution of mean size. Multiple comparisons were then conducted with the Hochberg method, for an adjusted p-value (Hochberg, 1988). An ANOVA was used to determine differences between sheltered and unsheltered behaviour (e.g., individuals in crevices or algae cover), based on individual sizes (normally distributed).

The second stage of the data analysis relied on empirical modelling to explain distribution patterns (density and mean size) according to environmental factors. To avoid misleading inferences due to similar effects of correlated variables (Morrissey and Ruxton, 2018), the collinearity between predictors to be included in the empirical models was assessed by applying a stepwise trait selection based on variance inflation factors (VIF), using the function vifstep() in the R package usdm (Naimi et al., 2014), with a VIF threshold set at th = 8, and checked with correlation plots. The selected predictors are described in Table 1. Apart from the variables excluded due to collinearity issues (dissolved solids, fine and medium sand fractions, and gravel), the two following variables were also left out. Salinity in the coastal area is more conservative than in estuarine environments. Since there is only one site at the estuary mouth (SE10) this could result in local salinity variations based on tidal regimes. Sea cucumbers are low mobility animals and salinity variations that are tidal dependent make a poor predictor for distribution patterns. Thus, other variables were chosen to represent a proximity to the estuary (e.g., longitude, turbidity, and ODO). Besides, longitude is not influenced by tidal regimes and salinity is, which, in turn, varies according to sampling hour and day (sampling logistics did not allow a standardised sampling regarding the tide). Shelter (presence or absence of crevices or seams where sea cucumbers can find shelter in), a binary variable representing a response to physical environmental stressors, like current, or predation, was also removed, because it was highly correlated to site current and caused convergence issues in the models. Sites without shelter were, simultaneously, the sites with the highest exposure to hydrodynamic influence (high current).

Table 1. Selected variables based on variance inflation factors to be included in empirical models assessing habitat preferences for Holothuria mammata.

The Generalized Linear Models (GLM) were implemented for two response variables: (1) density and (2) mean size, with Tweedie and Gamma distributions, respectively, and log-link functions (Tweedie, 1984). The ΔAIC was the measure derived from the AIC (Akaike information criterion) that was used as the model selection criterion and significant interactions were assessed considering the signs of the coefficients for the interaction and individual variables and with surface plots (Feld et al., 2016). A thorough description on the used models and the analytical procedures that supported the interpretation of the results are detailed in Supplementary Material 2. For relevant predictors in the models that showed a variations in time, a causality test, based on the ccf() cross-correlation function, was used to determine relationships between two time series: the relevant predictor variable and the response variable, across the sampling period. The autocorrelation function (ACF) plotted the results with all tested time-lags.

Results

Space-Time Distribution Patterns

The distribution range of H. mammata in the sampled area, like the other two commercial holothurians (Holothuria forskali and Holothuria arguinensis), was restricted to rocky reefs. Marine open areas without hard substrate had no occurrences and at rocky reefs no individuals were found beyond a maximum distance of 7 m from the transition border between rock and sand substrates. The estuarine area showed a wider heterogeneity concerning its environmental features and the dissimilarities from the adjacent marine environment were, mainly, related to the granulometric characteristics of the sediment, its organic content, turbidity, dissolved oxygen and temperature (Figure 2). The only estuarine site with H. mammata (or any of the aforementioned) was SE10, a site near the estuary mouth distinguished by the presence of a rocky outcrop and, in more than a year of sampling, no occurrence was observed at the two low depth nearby sites, SE9 and SE8, the latter in a seagrass meadow patch.

Figure 2. Principal Component Analysis (PCA) diagram of environmental data by sampled site, grouped according to habitat, classified as coast, estuary, and offshore.

Within the area of distribution of H. mammata 2114 individuals, ranging from 119 to 510 mm, were observed and counted in a total sampled area of 29,160 m². Density per replicate ranged from 0 to 120 ind/100 m² (or 12,000 ind/ha – a maximum value recorded at SC2) (Figure 3). The majority of individual counts were made on rock transects (94%), when compared to sand transects. Of those, 68% were settled on a sandy patch. In sand transects, the large majority was, in fact, settled on sand (82%) and remaining on small rocky patches. However, it should be considered that sand transects have very few rock outcrops and in rock transects sandy patches are common. Only a small number of individuals preferred a macroalgal cover, as only 5% occurred in the algal cover of rock transects and none in the algal cover of sand transects. Of all specimens, 34% were found to be sheltered in crevices at the time of sampling (daytime) and these individuals were overall smaller in size (ANOVA: p = 0.0086). When comparing sizes between substrates of settlement, individuals occurring on rock were also smaller than those on sand (ANOVA: p < 0.0001) and so were those occurring within an algal cover (ANOVA: p = 0.039).

Figure 3. Temporal density data for Holothuria mammata at Arrábida, Setúbal, during the nine month sampling period. Red dots at each month connected by a line represent the evolution of the mean.

Both density and mean size showed a significant spatial variation (KW: p = 2.474e-14 and p = 1.408e-13, respectively), but not temporal, although overall average density increased from spring up until early autumn (Figure 3). However, spatial differences in density were not related to site proximity (Figure 4). The most similar pair in density was SC5 and SE10 (KW, Hochberg: adjusted p = 1) and the most dissimilar pair was SC2 and SC3 (KW, Hochberg, adjusted p = 1.003e-10). Mean size, however, tended to increase towards the estuary, considering the four most distant sites (Figure 5). Size-classes presented a unimodal distribution at all sites. Mean size and density have a negative a linear relationship, with higher densities corresponding to lower size classes (R² = 0.12, F₁,₁₅₅ = 20.4, p = 1.24e-05).

Figure 4. Density of Holothuria mammata by sampling site, arranged according to distance to estuary (left to right). Red dotted line represents the overall mean value. * represents the mean of each sampling station.

Figure 5. Size-class distribution of Holothuria mammata by sampling site, arranged according to distance to estuary (bottom to top). Vertical lines represent the site mean.

Density Distribution Models

The GLM used to explain the density distribution of H. mammata, rendered only one final model with a ΔAIC > 2 and adjusted R² of 62% and an explained deviance of 68.6% (Table 2). Density, consistently higher at rock transects, decreased with increasing current, with pH (between the values of 7.71 and 8.54), with proximity to the estuary and increased with depth (Figure 6). The pH did not show a spatial variation, but a temporal one, increasing in the rainy months. The cross-correlation function showed a significant negative correlation at lag = 0 (months) (Figure 7A), with the peak of pH corresponding to the lowest average density value (Figure 7B). The significant interaction between depth and current, with a positive coefficient contrasting with the negative coefficients of its individual variables, showed an antagonistic effect, i.e., depth cancelled the negative effect of current on density. The interaction plot (Figure 8), detailing the effect, shows that between the depth values of 6 and 8 m, there was little variation in density, as opposed to lower depth values where the current had an effect on density, which decreased as current increased. In other words, the intensity of current only had a negative influence on density at low depths. At higher depths the hydrological characteristics of the site had less impact on density, as it even increased with increasing current at the two deepest sites.

Table 2. Fixed-effect GLM model results with a Tweedie distribution (R² = 62%), explaining density distribution patterns of Holothuria mammata at Arrábida, Setúbal.

Figure 6. Partial plots of the density model for all significant predictors – Fixed-effect GLM with a Tweedie distribution – explaining distribution patterns of Holothuria mammata at Arrábida, Setúbal.

Figure 7. Cross-correlation analysis for the time series of pH of the water column and mean density of Holothuria mammata, at Arrábida, Setúbal, during the nine month sampling period: (A) ACF plot depicting correlations (vertical lines) and significance threshold (horizontal dotted lines); (B) average monthly values for pH and density.

Figure 8. Surface plot depicting the interaction behaviour between current and depth in the density model, explaining distribution patterns of Holothuria mammata at Arrábida, Setúbal.

Size Distribution Models

The GLM used to explain the size distribution rendered four top models with a ΔAIC < 2. Model averaging for these models (full and conditional) revealed different levels of importance for several predictors, with substrate and distance to estuary (longitude) as the most relevant variables (Table 3 – full average), and coarse sand, ODO, TOM, Chla, Turbidity and the interaction depth:current as less important in explaining the variance of mean size for this species (Table 3 – conditional average). Considering the sum of weights of each variable in all four models, ODO and depth still presented a relatively relevant contribution in explaining the variance of the response variable. Coarse sand, TOM, the interaction depth:current and turbidity were relatively similar in contribution to the final model, but Chla was only represented in one model (Table 4). These results were corroborated by the assessment of the confidence intervals of the predictor’s coefficients (Figure 9). Nonetheless, of all the predictors regarded as less important, coarse sand and the interaction depth:current provided a more precise estimate of effect (narrower confidence intervals). Depth, although with a high importance (sum of weights), crosses zero and shows a wider confidence interval and, hence, not significant in either averaged model (full or conditional). Considering the interpretation of all model validation methods, the key predictors explaining size were substrate (larger individuals on sand) and longitude (larger individuals closer to the estuary). To a lesser extent, coarse sand (larger individuals associated to the higher granulometric fraction of sand), TOM (larger individuals at sites with more organic content in sediment: values between 0.3 and 8.5%), and turbidity (larger individuals at lower turbidity sites) also seem relevant. The relationship of larger individuals with more productive sites (Chla) was the least important association. The interaction depth:current was not assessed as to its type due to the unsignificant relevance of its individual predictors in the model. Of the predictors, ODO, TOM, and Chla, only the latter showed an obvious pattern, although like pH, it was not spatial, but temporal, with higher average values in the summer months. The cross-correlation function showed a significant positive correlation at lag = −1 (months) between Chla and mean size (Figure 10A), depicting an increase in local productivity preceded by a peak in mean length (Figure 10B).

Table 3. Fixed-effect GLM averaged model results with a Gamma distribution (averaged R² = 49%, as the average of all models weighted by their model weights), explaining size distribution patterns of Holothuria mammata at Arrábida, Setúbal.

Table 4. Model importance values extracted from the GLM averaged model, considering the four top models with a ΔAIC < 2.

Figure 9. Confidence intervals for the predictor coefficients of the model – Fixed-effect GLM averaged model with a Gamma distribution – explaining the size distribution patterns of Holothuria mammata at Arrábida, Setúbal.

Figure 10. Cross-correlation analysis for the time series of Chla and mean total length of Holothuria mammata, at Arrábida, Setúbal, during the nine month sampling period: (A) ACF plot depicting correlations (vertical lines) and significance threshold (horizontal dotted lines); (B) average monthly values for Chla and mean total length.

Discussion

The population of H. mammata at this NE-Atlantic rocky reef was the healthiest studied up to this moment. It showed the largest individuals reported so far, with a maximum total length of 51 cm, densities up to 120 ind/100 m² and unimodal size-class distributions. This reflects the good environmental conditions, successful recruitment and an unexploited population. The Mediterranean and other NE-Atlantic regions show populations of this species with lower maximum sizes, the majority between 25 and 38 cm (Navarro et al., 2013; González-Wangüemert et al., 2014, 2018; Marquet et al., 2017; Mustapha and Hattour, 2017; Siegenthaler et al., 2017). Previous records show that H. mammata in the NE-Atlantic attain bigger sizes than in the Mediterranean and the largest in the Portuguese coast, with Peniche, a highly productive marine area, presenting the previous largest documented specimen, with 43 cm (González-Wangüemert et al., 2016). Densities found in the present study were also the highest, compared to a maximum density of 55 ind/100 m² found in the Mediterranean, Aegean Sea (Aydin, 2019a), 16 ind/100 m² reported in the Canary Islands (Navarro et al., 2013), and particular low densities observed at Ria Formosa (south Portugal) of 1.2 ind/100 m² (Siegenthaler et al., 2017). However, densities tend to be much lower in habitats without rocky bottoms (Navarro et al., 2013) and Ria Formosa is a shallow mesotidal lagoon lacking the typical H. mammata habitat (Siegenthaler et al., 2017). Like at Peniche, the Arrábida coast is also a highly productive area, benefiting from upwelling events in the summer, nutrient inputs from the estuary and from a complex subtidal habitat (Wooster et al., 1976; Cabeçadas et al., 1999; Costa et al., 2013) that can have positive effects on the condition of the population. The high productivity together with the subtidal distribution of the species, unlike intertidal distributions (Siegenthaler et al., 2017), provide more stable environmental conditions, which may be the reason for higher abundances and lower variation in time.

This population was not only restricted to the marine environment but it preferred rocky reefs as opposed to other complex environments such as stands of macroalgae and seagrass meadows. Siegenthaler et al. (2017) had already reported that the preference of H. mammata for seagrass meadows at Ria Formosa did not reflect the species’ habitat preferences, but a distribution driven by the lack of rocky substrata and the consequent search for shelter and a more complex habitat (González-Wangüemert et al., 2014, 2016). Despite the higher occurrence found on rocky substrate, smaller animals were the ones sheltered in crevices or macroalgae, contrary to those on sand. There are a number of sea cucumber predators, with sea stars as the most common, particularly at the juvenile stage (Francour, 1997). Hence, those larger that venture more into open spaces (for feeding activities) are less prone to be preyed upon, as predation risk declines with growth (Shiell and Knott, 2008; Purcell, 2010; Ceccarelli et al., 2018). It is important to stress that while sandy patches in rocky areas have shelter close by, sand areas do not, which may be the cause for the proportion of individuals on sand to be higher when considering the individual presence, rather than the transects density (sandy areas vs. rock). So, this may be either a defence strategy for smaller individuals, that use sandy patches to feed, and an opportunistic behaviour exclusive of larger animals that move more into sandy areas, or instead, a diel feeding behaviour of lower size-classes that may be more active at night, as this is a species with reported nocturnal activity (Borrero-Pérez et al., 2010; Navarro et al., 2013; Siegenthaler et al., 2017). Nighttime was not covered in the present study, but there was a large proportion of active individuals at this NE-Atlantic rocky reef during the day. This was not only assessed by the number of individuals on sand, because H. mammata feeds exclusively on sand in this area (Azevedo e Silva et al., 2021) and needs to move from the rock substrate to a sandy area or patch to feed, but also demonstrated by the observation of feeding behaviour at sampling time for several individuals (footage obtained at the study area between 10 and 12 h a.m. at https://seacucumber.eu/en/o-pepino-do-mar).

Like other sea cucumber species (Džeroski and Drumm, 2003; Mendes et al., 2006; Eckert, 2007; Shears and Babcock, 2007; Dissanayake and Stefansson, 2010; Domínguez-Godino and González-Wangüemert, 2020), the population studied in the present work showed a patchy distribution, with areas of high abundance followed by zero or near zero sites. However, this patchiness is size dependent. Lower size-classes aggregated in higher density groups in more stable environments (away from the estuary).

The models results suggest that abiotic and biotic parameters of the water column are not the main drivers shaping the distribution of H. mammata in the area. Similarly to other sea cucumber species (e.g., Sloan and von Bodungen, 1980; Sonnenholzner, 2003; Mendes et al., 2006; Morgan, 2011; Domínguez-Godino and González-Wangüemert, 2020), its habitat preferences depend on environmental stability, as opposed to the estuarine variability, low hydrodynamic conditions and the presence of rocky substrate that may offer shelter and habitat complexity, more likely to increase food quality. Depth plays an important role by providing a quieter refuge in areas that are more susceptible to tidal influence. Individuals tend to use bathymetry to their advantage in areas where current is stronger, occupying deeper areas. The significant physical-chemical predictor of the water column explaining density was pH. However, pH shows little variation in space, but is highly variable in time. Its increase in the autumn and winter months is most likely related to continental run-off from the limestone ridge of Arrábida that characterises the entire coastline (Kullberg et al., 2012). This temporal variation showed a negative correlation with the mean density that, within this time-frame, represents an aggregation phenomenon that increases capturability during the months that coincide with the pH decrease in this coastal area. In echinoderms, pH may affect abundance, distribution, and also reproductive success, as a consequence of induced physiological stress (Lawrence and Herrera, 2000; Walag and Canencia, 2016). Coincidently, the highest densities found occurred during spring/summer months, which correspond to the species reproductive season (unpublished data for the area; Marquet et al., 2017), and may reflect a relationship between reproduction and suitable environmental conditions, where pH can play a role, promoting an aggregation behaviour during this period (Hamel and Mercier, 1996b; Mercier and Hamel, 2009; Leite-Castro et al., 2016). The model results corroborate the relationship between density and size. Despite the lower densities on sand, these were represented by the largest individuals. Smaller individuals were distributed in rocky substrate and at an increasing distance from the estuary, where, by comparison, there is less environmental variability. The relationships with H. mammata sizes and dissolved oxygen, organic content in the sediment and Chla were, in general, difficult to assess due to the low importance of these predictors in the model as well as their inconspicuous patterns. Only Chla showed an evident temporal variation, increasing in the summer, which is a common phenomenon in the area, due to seasonal upwelling cycles (Wooster et al., 1976). Similarly to pH and density, this rising local productivity is associated an aggregation behaviour of larger size-classes in the spring/summer months. Hence, there is a strong suggestion of an aggregation pattern of breeders in a specific period, which relates to environmental conditions that potentially favour the physiological process of gonad development. The species’ preference for coarser sand fractions may be related to a more diverse diet and less related to organic sediment content (Underwood and Barnett, 2006). So, these associations may be explained by the displacement ability of larger individuals towards more favourable conditions, here expressed by coarser sediment fractions, higher productivity, less turbidity and more organic content of the sediment (never higher than 8.5%), at a lower risk of predation or due to a higher resilience, both privileged by size.

Conclusion

In general, H. mammata prefers a more stable environment, which is more important for smaller individuals, while bigger size classes tend to move to less stable environments, probably in an opportunistic fashion due to their size related resilience. This behaviour could be related to feeding conditions associated to lower predation susceptibility of bigger animals, or with a reproductive behaviour leading the broodstock, as they reach maturity size, towards more favourable areas for reproduction, since sea cucumber larvae are pelagic (Venâncio et al., 2021) and require currents for dispersion (Pedrotti and Fenaux, 1992).

Conducting research in an unexploited area provided important advantages in the study of biological and ecological traits, closer to a pristine condition, in comparison to populations that are subject to harvest, several of which already showing signs of exploitation by a reduction of larger size classes, thus creating multimodal distributions (e.g., González-Wangüemert et al., 2014, 2018; Aydin and Erkan, 2015). Although, this work was conducted in a limited area of an Atlantic rocky-reef, there is a consistency in the drivers that shape the distribution of different detritivorous sea cucumbers, when compared to other studies (e.g., Sloan and von Bodungen, 1980; Sonnenholzner, 2003; Mendes et al., 2006; Morgan, 2011; Domínguez-Godino and González-Wangüemert, 2020). These species tend to prefer the available shelter, a more complex habitat for feeding and areas with less environmental variability. This provides guidance for future stock management measures in a time when only a small proportion of the species’ distribution range is studied. At the studied rocky coastal area, the habitat preferences of H. mammata showed that, at this point, even lacking important information like growth parameters or size at first maturity to allow the establishment of recovery periods for the population and determine, e.g., fishing quotas for each size-class, some red flags can be raised as to the potential susceptibility of H. mammata to fisheries. The patchy distribution and the high density areas should be regarded as sensitive areas, as these are more attractive for fishermen, allowing a reduced fishing effort and increased yields, which can lead to a significant reduction in spawning biomass. Access to high density of breeders should also be restricted to allow a spillover effect and help sustain the adjacent populations (Purcell and Kirby, 2006). On an undesirable note, and to the best of our knowledge, no study, so far, covers the distribution of wild juveniles of H. mammata under 100 mm. The understanding of juvenile habitat preferences, distribution and the integration of these size-classes in growth models is key knowledge and, currently, a large gap. During the sampling campaigns no small juveniles were observed, which suggests a different habitat preference or, most likely, a cryptic or diel behaviour, as observed for other sea cucumber species (e.g., Sloan, 1979; Shiell, 2004; Purcell, 2010; Soliman et al., 2019).

Although, currently, the exploitation of NE-Atlantic sea cucumbers is restricted in space and to a relatively low number of illegal fishermen, the high-income potential that these animals yield will soon widen the interest that easily leads to a high fishing pressure. This study represents a fundamental step towards a better understanding of the populations’ dynamics and can contribute with vital knowledge support decision-making processes on stock management.

Data Availability Statement

The datasets presented in this article are not readily available because these are under an embargo period. Requests to access the datasets should be directed to PF, pmfelix@fc.ul.pt.

Author Contributions

PF, AB, and JC designed the study. FA, TS, and PF carried out the field work. PF, TM, and CR developed the data analysis. AP, RM, JS, and EV provided the contributions to several aspects of the study and the manuscript. All co-authors revised the manuscript.

Funding

This study was financed by the Operational Programme Mar2020, MAR-02.01.01-FEAMP-0052, “Newcumber – Avanços para o cultivo sustentável de pepinos do mar”. It received further financial support from Fundação para a Ciência e Tecnologia (projects UIDB/04292/2020, UIDB/00006/2020, CoastNet – PINFRA/22128/2016, AB with the Scientific Stimulus Programme – CEECIND/00095/2017 and FA with the individual research grant 2020.09563.BD). This publication was financed by Portuguese national funds through FCT – Fundação IP under project reference UIDB/04292/2020, and by the European Union’s Horizon 2020 Research and Innovation Program under grant agreement N810139: Project Portugal Twinning for Innovation and Excellence in Marine Science and Earth Observation – PORTWIMS.

Conflict of Interest

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

Publisher’s Note

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

Acknowledgments

The authors would like to thank Ricardo Lima for developing the map for the study area and Pedro Segurado for his crucial help with the interaction plots.

Supplementary Material

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

References