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

There is an increasing demand for sea cucumbers, for human consumption, mainly from Asian markets and, as a consequence, NE-Atlantic species are now new targets for exploitation and exportation. Holothuria mammata is one of the most valuable species in Europe. However, the lack of historical economic interest in this species in most European countries has also led to a lack of studies concerning biological and ecological aspects on wild populations and this is a major issue for stock management. 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, as a function of environmental conditions. For that, a population from a marine protected area was followed for 1 year at 1.5-month intervals. Throughout the coastal area, six sites were selected and at each sampling campaign three random transects per site and substrate (rock and sand) in which all H. mammata individuals were counted and measured. For each site and survey several environmental parameters of interest, from the water column, the sediment and substrate cover, were also measured. Generalized Linear Models were used to model the spatial and temporal distribution of the species according to environmental conditions, to determine the species’ habitat preferences. The distribution models indicate that abiotic and biotic parameters of the water column are not the main drivers shaping the distribution of H. mammata. The species has a patchy distribution, and its habitat preferences depend on environmental stability, the presence of shelter and habitat complexity, which is more important for smaller, more vulnerable, individuals, while bigger size classes tend to venture more into less stable environments in an opportunistic fashion. The knowledge of these population traits is determinant to develop stock management measures, which are now urgent to prevent the depletion of commercial sea cucumber populations in the NE-Atlantic. Sustainable fisheries policies should be developed and start by considering to delimit fishing areas and periods, considering the species spatial and temporal distribution 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 physicochemical 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., 2011Schneider et al., , 2013Purcell 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 sizeclasses, 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., 2015Siegenthaler et al., , 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., 2012Purcell et al., , 2016Kashio 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., 2015Siegenthaler et al., , 2017Marquet 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., 2014Dereli 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., 2014Aydin 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, 2019aAydin, ,b, 2020Sousa 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.

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.

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 FIGURE 1 | Map of the study area, with sampling sites, grouped according to habitat type, classified as; estuary (yellow); coastal (red); and offshore (blue).
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 Frontiers in Marine Science | www.frontiersin.org Longitude was chosen as a proxy of environmental stress, as it represents distance to estuary mouth (longitude increases with distance). Estuarine environments have more environmental variability than its adjacent coastal areas.
Current Continuous m.s −1 Current represents the hydrodynamism of each site, as an environmental stressor. Reflects the tolerance of sea cucumbers to hydrodynamic conditions when related to density or size class.
Coarse sand Continuous % A granulometric feature of the sediment that may be related to feeding preferences. Of all sand fractions, the coarse sand was that which best explained the granulometry of the sediment, by correlating with all others.
Silt Continuous % A granulometric feature of the sediment that may be related to feeding preferences. Silt represents the finer fraction (<63 µm) of the sediment.
Total organic matter (TOM) in the sediment Continuous % Total organic content of the sediment is a measure of organic content availability that may be related to feeding preferences, although not providing information on nutritional quality.  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). 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.

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.
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 2 . Density per replicate ranged from 0 to 120 ind/10 m 2 (or 12,000 ind/haa 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).
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 2 = 0.12, F 1 , 155 = 20.4, p = 1.24e-05).

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 2 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 crosscorrelation 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.

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

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/10 m 2 and unimodal size-class distributions. This reflects the good environmental conditions, successful recruitment and an unexploited population. The p-values significance codes: 0 "***" 0.001; "**" 0.01; "*" 0.05; "." 0.1; " " 1.
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., 2014Marquet 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/10 m 2 found in the Mediterranean, Aegean Sea (Aydin, 2019a), 16 ind/10 m 2 reported in the Canary Islands (Navarro et al., 2013), and particular low densities observed at Ria Formosa (south Portugal) of 1.2 ind/10 m 2 (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 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. 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. 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 sizeclasses that may be more active at night, as this is a species  p-values significance codes: 0 "***" 0.001; "**" 0.01; "*" 0.05; "." 0.1; " " 1.
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  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., 2014Aydin 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 sizeclasses 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.