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.

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.

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.

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.

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.

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.

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.