Introduction

Sea cucumbers play an important role in marine ecosystems (^{Purcell et al., 2016}), and present outstanding commercial value throughout the world (^{Purcell et al., 2018}), however little is known about their populations in Brazil (^{Mendes et al., 2006}; ^{Ponte & Feitosa, 2019}; ^{Souza Jr. et al., 2017}; ^{Ventura et al., 2013}). The greed for these high priced-organisms has led to intense capture in more than 70 Countries, threatening the populations of various species around the planet (^{Purcell et al., 2013}; ^{Purcell et al., 2023}), including Latin America (^{Conand, 2018}; ^{Sonnenholzner, 2021}).

Although a wide diversity of sea cucumbers is described for shallow waters of the Brazilian coast (^{Martins, 2012}; ^{Martins & Souto, 2020}; ^{Martins & Tavares, 2021}; ^{Martins et al., 2022}), only three species are described for southern Brazil (^{Rupp et al., 2023}). Among them, Holothuria Halodeima grisea Selenka, 1867 is considered the most common, being recorded in the intertidal region at the base of the rocks, usually in contact with the bottom sand (^{Bueno et al., 2015}; ^{Martins, 2012}; ^{Mendes et al., 2006}; ^{Tiago & Ditadi, 2001}; ^{Tommasi, 1969}). Nevertheless this species is already threatened by unreported and unregulated capture in Brazil (^{Ponte & Feitosa, 2019}; ^{Rupp & Marenzi, 2021}; ^{Souza Jr. et al., 2017}). The geographic distribution of H. (H.) grisea is predominantly tropical, and includes the west coast of Africa, Florida, Gulf of Mexico, Panamá, Caribbean islands, Colombia, Venezuela, and Brazil up to Santa Catarina (^{Martins, 2012}; ^{Pawson et al., 2010}). In latter location ^{Rupp et al. (2023}) recorded the occurrence of H. (H.) grisea in nine out of the eleven intertidal sampled sites, with higher densities (> 1.7 ind.m^-2) recorded in the central portion of the state, and none found further south than the locality of Garopaba.

The objective of the present study was to evaluate the abundance and densities (individuals.m^-2) of H. (H.) grisea within a spatial-temporal perspective, at a locality where a previous study has been conducted (^{Mendes et al., 2006}) nearly 16 years earlier. Additionally, biometric parameters of the population as such length and weight distributions were determined, as well as growth parameters. Such basic information could be useful for future conservation practices as well as for fishery management and aquaculture of this species.

Materials and methods

Field methods: The study site is in the littoral zone of Armação Itapocoroy cove, Penha, Santa Catarina (26º 47’S; 48º36’W) (Fig. 1), which is a low-energy wave action area. The site has a slight slope and rocky-sandy substrate consisting of coarse grain-sized sediments (^{Mendes et al., 2006}). This site is located approximately 150 km from the southernmost point where the Holothuria (H.) grisea have been recorded (^{Rupp et al., 2023}). Four Holothuria (H.) grisea sampling campaigns were carried out between August 2019 and October 2020 during low spring tides. Restrictions imposed by the Covid-19 pandemic prevented sampling from being carried out from autumn to winter 2020. On each occasion, three two-meter wide transects (T1, T2 and T3) distant 15 m from each other (Fig. 2) were deployed perpendicularly to the shore, with lengths ranging from 33 to 50 m (Table 1), depending on the tidal level. The transect areas were surveyed by visual and tactile inspections and the numbers of H. (H.) grisea in each one-meter division of the transects were recorded. The samplings always started at T2, then proceeding to T3 and finally to T1. In three opportunities it was not possible to finalize the sampling at T1 and T3 due to tidal surge. The rocks at the sampling area were not revolved in order to avoid habitat disturbance. The length of the transects was divided in three sectors based on tidal exposure: A - upper intertidal (length 0 to 16 m), B - lower intertidal (length 17 to 33 m) and C - subtidal (length 34 to 50 m). Their depths were determined in the middle of each sector. To do that, tidal tables from the Brazil’s Navy Hydrography and Navigation Center (^{Centro de Hidrografia e Navegação (CHN), n.d.}) and from a nearby tide gauge from the Santa Catarina State Institution for Agricultural Research and Rural Extension (Empresa de Pesquisa Agropecuária e Extensão Rural de Santa Catarina, n.d.) were used as references. The rugosity of each sector was estimated using the chain method (^{Luckhurst & Luckhurst, 1978}). A metal chain with small links was laid along the central portion of the transects, so that it followed the contours of the substratum as closely as possible. The ratio between the contour length of the chain and the sector actual length (Rg) was used to infer rugosity index:

Ri = 1-(1/Rg)

Where Ri is the rugosity index and Rg is the ratio between the chain length and the actual length of the sector of the transect (^{Mendes et al., 2006}).

Fig. 1 Location of the study site at Armação do Itapocoroy, Penha (SC). 

Fig. 2 Plots of the transects at Armação do Itapocoroy, Penha (SC) (A - Google Earth); (B) aerial photo of the sampling area. (Arrows indicate North). 

A group of 90 animals was randomly sampled, on each campaign, for in loco biometry. The following parameters were registered: total length and width with contracted body, and total weight after light pressure to expel the liquid from the internal cavity. Linear measurements were carried out with a precise ruler (± 0.5 mm) and a field DC digital scale was used to determine weights (± 0.1 g). Samplings were non-destructive and the animals were immediately returned to their habitat after biometry. The moon phase and the following environmental parameters were recorded in each sampling occasion: air temperature, sky coverage (clouds), water temperature and salinity. Tidal information was obtained at CHN tables (CHN, n.d.). This study was carried out under license for collection of native fauna of the Brazilian Ministry of Environment MMA/SISBIO No. 68215.

Statistic approach: Normality and homoscedasticity of the datasets were verified respectively with Shapiro-Wilk and Levene’s tests using Statistica^R, version 12 (Statsoft Inc. USA), and density data was square root transformed in order to meet those premises (^{Zar, 1999}). Two-Way ANOVA was used to compare densities on spatial-temporal scale (among transects, sectors and season) (^{Sokal & Rohlf, 2000}). Significant differences were considered when P < 0.05 and Tuckey’s test was used as post-hoc analysis. The comparison of mean length and weight of Holothuria (H.) grisea among the sampling campaigns was explored using One-Way ANOVA, however the assumption for this analysis were not met, even after data transformation, therefore the non-parametric Kruskal-Wallis test was used in lieu of ANOVA (^{Sokal & Rohlf, 2000}).

Considering that it was not possible to complete Sector C on T1 in 3 opportunities and on T3 in one opportunity, only the sectors fully sampled were included in the statistical analyses (Table 1). The comparison of the mean densities on temporal scale was carried out considering T2 (sectors A, B and C) for all sampling occasions, as well as for all transects using sector B only. The comparison among all sectors and transects was carried out within the summer-20 campaign and the comparison of the densities among sectors B was carried out considering all transects and seasons.

Aiming to explain the density of sea-cucumbers based on both the mean depth of the sectors and their rugosity indexes, linear regression analyses were carried out. The density of animals was the response variable, while the depth and rugosity were tested as explanatory variables, both separately and combined as a multiple regression model, using the software “R”, version 4.3.0 (^{R Core Team, 2023}). The relation between size and weight of the studied population was explored using the power model regression:

Y = aX^b

Where Y = weight, X = length, b = allometric coefficient, a = intercept (^{Gould, 1966}).

Considering that body measurements in holothuroids is often imprecise due to body wall elasticity, the compound index which combines body length and width (SLW = square root of the length-width product) was used to further analyze its relationship with body weight, as indicated for several species of sea-cucumbers (see ^{Poot-Salazar et al., 2014} for details). To increase growth parameter accuracy, lengths were recalculated using a power regression between length and SLW (a = 1.0731, b = 0.6855).

Growth analyses of different SLW corrected length classes (Le) were conducted by Electronic Length Frequency Analysis (ELEFAN) in R 4.3.0 (^{R Core Team, 2023}) using the TropFishR package (^{Mildenberger et al., 2017}). The von Bertalanffy growth function was applied as follows:

Where Le_t is the predicted corrected size at age t; Le_∞ is the asymptotic size; k is curvature parameter per year, expressing the rate at which Le_∞ is approached; and t_o is the theoretical ‘‘age’’ if the organism were to have a size equal to zero. Relative ages at certain lengths and weights were generated by replacing the growth parameters with their respective values as follows:

Results

Environmental variables: The environmental variables recorded during the surveys carried out from August 2019 to October 2020 are shown in Table 2. The sectors’ surveyed areas, and the respective depths and rugosity indexes are presented in Table 3.

Densities and spatial-temporal distribution: The total number of Holothuria (H.) grisea counted during the study was 566 individuals out of which 360 were used for biometrics. The highest density of H. (H.) grisea recorded per unit area of the transects was 22.5 individuals*m^-2 (T3-winter-09). The densities varied along the transects for all sampling campaigns (Fig. 3).

Fig. 3 Densities of Holothuria (H.) grisea (ind*m^-2) along the transects for the sampling campaigns carried out in winter-2019, spring-2019, summer-2020 and spring-2020. Uppercase letters indicated the transect sectors: A. upper intertidal, B. lower intertidal and C. subtidal. 

The mean densities of Holothuria (H.) grisea within the different sectors of each transect for all sampling campaigns are presented in Fig. 4. It clearly shows that transect A (upper intertidal) had lower densities than sectors B (lower intertidal) and C (subtidal) in all transects and samplings. The mean densities for the samplings in Winter/2019, Spring/2019, Summer/2020, and Spring/2020 considering all the surveyed sectors and transects were 2.11, 1.21, 1.67, 1.94 individuals*m^-2, respectively. However, when only the lower intertidal and subtidal sectors are taken into account (excluded sector A where specimens were rarely found), the mean densities for winter/2019, spring/2019, summer/2020, and spring/2020 were 3.35, 2.02, 2.38 and 2.80 individuals*m^-2, respectively.

Fig. 4 Densities of Holothuria (H.) grisea for each sector for all the transects and sampling campaigns. The boxes indicate the interquartile range divided by the median and the whiskers indicate the lowest and highest results. A. upper intertidal, B. lower intertidal and C. subtidal. 

The comparison of the densities among transects, sectors, and seasons are presented in Table 4 (Two-Way ANOVA). The analyses indicate no difference among the surveyed seasons, in any of the compared sectors and transects. On the other hand, significant difference among sectors was recorded in transect 2 (considering all seasons), with significantly higher densities recorded in sector C, intermediate in sector B, and lowest densities in sector A (Tuckey test, p > 0.001). In Summer-19 densities recorded on sector C were also significantly higher than sector B and A, respectively (Tuckey test, p < 0.001), but there were no differences among transects. Considering sector B only, when all transects and seasons were considered, there was significant difference among transects, with lower densities recorded in transect 3 (Tuckey test, p < 0.01).

The effect of bathymetry and rugosity on the density of animals: Regression analyses evidenced that the density of sea-cucumbers is significantly correlated with both the depth (p <0.001) and the rugosity (p <0.001) of the studied sectors (Fig. 5). The levels of explained variance, however, differed significantly, being 75 % and 35 %, respectively for models considering depth (Residual standard error: 0.7 ind.m^-2) and rugosity (Residual standard error: 1.18 ind.m²) as explanatory variables. The former and best model evidence that the density of sea cucumbers tend to increase by 1.3 individuals.m^-2 for the increment of every 10 cm in depth (density = -5.66+13.5*depth). A gain of 4 % in terms of explained variance is obtained for the multiple linear model considering both depth and rugosity as explanatory variables (R² = 0.79, p <0.001, model: Density = -5.4+ 11.88*Depth + 34.16*Rugosity).

Fig. 5 Regression lines considering density of Holothuria (H.) grisea as response variable and A. rugosity and B. depth as explanatory variables. 

Length and weight distributions: The length of the measured Holothuria grisea specimens ranged from 5.2 to 22.5 cm and their weight ranged from 6.0 to 230 g (Fig. 6). The mean length of the population when data from all sampling occasions were combined was 12.54 cm (SD = 2.53) (mode = 13 cm; median = 12.4 cm). Approximately 75 % of the sampled specimens were between 10 - 14 cm. The overall mean weight was 60 g and nearly 60 % of the population was between 55 and 93 g. The relationship between body length and weight of H. (H.) grisea (Fig. 6) is well represented by the power model (p < 0.001) with the variation in length explaining 62 % of the weight variability. When the same database is analyzed separately by different sampling occasion, no difference in length (Kruskal-Wallis, H = 5.27, p = 0.152) and weight (Kruskal-Wallis, H = 3.31, p = 0.345) was observed (Fig. 7). Therefore, growth analyses were represented for the entire dataset.

Fig. 6 Length and weight histograms of Holothuria (H.) grisea as well as the length-weight regression pooled for all the sampling campaigns. 

Fig. 7 Length and weight histograms of Holothuria (H.) grisea as well as the length-weight scatterplots recorded on the surveys carried out in winter 2019 (WI-19), spring 2019 (SP-19), summer 2020 (SU-20) and spring 2020 (SP-20). 

Growth analysis: The relationship between the compound index SLW and body weight of Holothuria (H.) grisea is presented in Fig. 8A. The SLW calculated for Holothuria grisea specimens ranged from 2.79 cm to 9.84 cm. The SLW calculated explained 75 % of variation in body weight, and was a more effective predictor of weight than length (62 %).

Fig. 8 A. Relationship between the compound index SLW (square root of the length-width product) and body weight (g) of Holothuria (H.) grisea. B. Individual growth of Holothuria (H.) grisea expressed by the von Bertalanffy growth function in cm*y^-1 fitted to the corrected lengths (Le), and the relative ages generated by replacing the growth parameters calculated by ELEFAN in the inverted von Bertalanffy growth function. 

The ELEFAN method (Electronic Length Frequency Analysis) yielded the best-estimated growth parameter results: Le_∞ = 13.63, k = 0.72 y^-1, t_o = 0.79. Graphic representations of the growth curve of Holothuria (H.) grisea showed that the youngest sea cucumber was one year-old, and the oldest had approximately two years and a half of age. Growth rates tended to decrease between the first and the second year of life (Fig. 8B).

Discussion

The results of the present study evidence that the densities of Holothuria (Halodeima) grisea observed on the different sectors of the transects are positively correlated with depth, and therefore with the period and frequency at which these areas are subject to air exposure due to the tidal level variation. The coast of Santa Catarina is characterized by a mixed microtidal regime with predominance of astronomic semi-diurnal tides, which means that there is a complete tidal cycle every 12 hours, and a complete lunar cycle every 29.5 days (^{Salles-de Araujo, 2020}). Tides are also frequently influenced by meteorological factors which tend amplify the predicted astronomical tidal range (^{de Souza & Schettini, 2014}; ^{Salles-de Araujo, 2020}). During new and full moons, the tidal amplitudes are nearly 1.2 m (spring tides), whereas during crescent and waning moons the neap tides display low amplitudes (~ 0.5 m) (^{Salles-de Araujo, 2020}). In this manner, the sea cucumbers in the upper intertidal level are subject to frequent air exposure during the tidal cycles, whereas in the subtidal sector, the frequency of air exposure is negligible, and in the lower intertidal sector the periods of air exposure tend to occur only for a few hours during spring tide periods (full and new moons). Air exposure and desiccation represent a stressful situation for sea cucumbers (^{Hou et al., 2019}) and therefore they would tend to avoid the areas prone to long periods of emersion.

The similarity in densities observed among the transects when all sectors were considered is explained by the overall similarity in the depth range and rugosity indexes when considering the total area of the transects. On the contrary, comparing sector B alone, the lowest density recorded in the transect 3, is clearly related to the lower rugosity of this segment, since the depth of sector B among transects is similar. Therefore, within similar depths, the rugosity tends to be an important factor affecting the aggregation of H. (H.) grisea. In Canary Islands, ^{Navarro et al. (2013}) also found significant correlations relation between density of sea cucumbers and rugosity of the substrate.

^{Mendes et al. (2006}) found that H. (H.) grisea occurred more densely in the subtidal stratum and displayed an aggregated distribution pattern, which was correlated to areas with high substrate rugosity, which on its turn, was highly correlated with rock coverage. Such observations were corroborated by the present study, and furthermore, we demonstrate that depth is a primary factor explaining the observed density patterns while rugosity of the substrate is secondary. ^{Leite-Castro et al. (2016}) observed a clear correlation between aggregative behavior and the reproductive cycle of H. (H.) grisea in northeast Brazil, with strong negative correlation between aggregation and gonad maturity. The highest degree of aggregation (groups of > 3 individuals under the same rock) occurred when the gonads were poorly developed or post-spawned. Conversely, the highest frequency of solitary individuals occurred when most of them were mature. In that region, the reproductive peak of gametogenic activity is from December to February when aggregative behavior was minimal. In southern Brazil, information about the reproductive cycle of this species is still limited. The only study available (^{Bueno et al., 2015}) suggests that mature individuals are found year round with increased sexual activity in February, however this study was carried out with limited samplings and organisms. The reproductive pattern of H. (H.) grisea, may affect its aggregation behavior in southern Brazil, however further investigations on reproduction of this species are required to examine such hypothesis.

Information about the status of the wild populations of Holothuria (H.) grisea is scarce. ^{Souza Jr. et al. (2017}) refers to densities up to 40 ind.m^-2 in the locality of Bitupitá (Ceará) in northeast Brazil before the onset of sea cucumber fisheries, whereas subsequently to the establishment of fishing pressure, ^{Farias-Dias (2012}) records mean densities of 0.54 ind.m^-2 in the same area, which indicates a dramatic drop (nearly 98 %) in sea cucumber densities. The characterization of that fishery is presented by ^{Ponte and Feitosa (2019}), which indicates that in some areas the catch could reach up to 6 600 Kg*yr^-1. However, no information is presented about the status of the wild populations on those fishing sites. ^{Rupp et al. (2023}) found densities ranging from 0.11 to 3.66 ind.m^-2 at the intertidal region on different sites along Santa Catarina coast, and Rupp and Marenzi (2021) inform the occurrence of eventual unregulated and unreported fishery of sea cucumbers also in the state of Santa Catarina. ^{Mendes et al. (2006}) reports mean densities ranging from 0.74 to 8.04 ind*m^-2 in the same region of the present study. Although that study does not present the absolute densities recorded per sampled metric unit as in the present one (Fig. 3), a comparison of the mean densities (ind*m^-2) presented in that study, carried out 16 years before the present one, is interesting. The mean densities obtained in the former study at the lower intertidal and subtidal strata were respectively 3.05 and 6.65 ind*m^-2; and for winter-03, spring-03 and summer-04 the mean densities were respectively 5.68; 5.34 and 4.14 ind.m^-2, respectively. In the present study the mean densities recorded in the lower intertidal and subtidal regions were 2.17 and 3.31 ind*m^-2, respectively; and for winter-19, Spring-19 and summer-20 were 3.35, 2.02, 2.38 ind*m^-2, respectively. The comparison of both studies suggests that the densities recorded in the present study are remarkably lower than those recorded in the previous one for equivalent strata and seasons (28 % to 62 % drop depending on strata end season). However, it is not possible to infer whether the differences from both studies are due to biological factors or due to an eventual effect of undercovered fishing activity. In any case, differently from the northeastern Brazil, where a dramatic drop in the population densities were found after the onset of the sea cucumber fisheries (^{Souza Jr. et al., 2017}), the reduction in densities of H. (H.) grisea observed in Santa Catarina does not seem to be so striking, but it certainly rises a concern about eventual anthropic influence on wild sea cucumber stocks.

Juveniles of most sea cucumber species are seldom observed in great numbers in the wild (^{Shiell, 2004}; ^{Shiell, 2005}) as also noted for H. (H.) grisea by ^{Leite-Castro et al. (2016}) in northeast Brazil. In the present study there was a low frequency of juveniles smaller than 8 cm and no individual was smaller than 5.2 cm. The reasons for such pattern in holothuroids comprise two main hypotheses: that young juveniles occur alongside with adults but displaying a cryptic behavior, therefore being difficult to find (^{Cameron & Fankboner, 1989}; ^{Rogers et al., 2021}); or else, juveniles settle in specific habitats outside adult areas (^{Mercier et al., 2000}a; ^{Mercier et al., 2000}b). In fact, ^{Aquino-Souza and Gomes-Filho (2023}) found high densities of H. (H.) grisea juveniles (< 7.5 cm) associated with seagrass beds (Halodule wrightii) in northeast Brazil, favoring the latter hypothesis, however only one juvenile below 5 cm was found. Consequently, the habitat for juveniles smaller than 5 cm is still to be found for H. (H.) grisea along its distribution range. Whereas seagrasses forms extensive meadows in northeast Brazil, they are rare in Santa Catarina (^{Marques & Creed, 2008}) and Halodule wrightii is non-existing, therefore other habitats should be searched as settlement area for H. (H.) grisea. However, it must be pointed out that, in the preset study, we avoided substrate disturbance, and rocks were not removed and sediments not revolved, therefore the hypothesis that the juveniles are more cryptic, hiding deeper under the rocks, crevices and cracks, could not be ruled out at this moment.

The temperatures recorded in the present study (both air and water) are in accordance with the expected seasonal variation in this region. Salinity displayed little variation among the campaigns, as all of the samplings were carried out with fair weather, even though severe rainstorms may often occur, decreasing the salinity in coastal areas of Santa Catarina (^{de Souza et al., 2016}). All the surveys of the present study were carried out in periods of full and new moons, when the lower tidal levels are recorded. Although studies demonstrate that moon phase is an important factor affecting holothurian reproduction (^{Mercier et al., 2007}), it is not clear whether it directly influences sea cucumber abundance and densities.

Sea cucumber fishery is an important activity in northeast Brazil generating income for fisherman of coastal communities, however it is an unregulated, not declared, and unmanaged activity, and a severe decrease of wild stocks have been reported (^{Leite-Castro, 2016}; ^{Ponte & Feitosa, 2019}; ^{Souza Jr. et al., 2017}). ^{Ponte and Feitosa (2019}) provide information about the size distribution of the catch in certain localities in the state of Ceará. The average length of captured H. (H.) grisea varied among fishing sites, ranging from 10.93 cm (SD = 2.01, median = 10.80 cm) in the Conflict Zone to 14.43 cm (SD = 2.36, median = 15.67 cm) in Caraúbas, with and average size among sites of 12.25 cm. Additionally, ^{Farias-Dias (2012}) reports a mean body length of 12.59 cm (SD = 2.42) in Bitupitá. The mean length of H. (H.) grisea found at Armação do Itapocoroy in the present study (12.54 cm, SD = 2.53, median = 12.4 cm) suits within the range of those reported for northeast Brazil. It was not possible, however, to compare the weights of the sea cucumbers from both regions, since the information presented by ^{Ponte and Feitosa (2019}) refers to the eviscerated organisms, whereas our study evaluated the whole body weight.

The relationship between body length and weight of H. (H.) grisea is well represented by the power model (P < 0.000001) with the variation in length explaining 62 % of the weight variability, and the obtained “b” value was 1.8765. ^{Farias-Dias (2012}), using the same model found a determination coefficient of 54 % and a “b” value of 1.225 for the same species in northeast Brazil. The obtained “b” value, different from 3, indicate negative allometry (^{Gould, 1966}), which means that the body length of the species grows faster than its weight. Negative allometry is commonly observed in other sea cucumbers (^{Khodja & Mezali, 2023}; ^{Poot-Salazar et al., 2014}). The regression of the compound index SLW (square root of the length-width product) and weight indicates that the use of this index enhance the weight estimation accuracy by 75 % rather than the Le. On the other hand, the determination of Le allowed the calculation of growth parameters and age determinations. The curvature parameter per year (K) also approximate other determinations conducted for sea cucumbers in Mexico and Portugal (0.20-0.88) using ELEFAN (^{Olaya-Restrepo et al., 2017}; Poot-Salazar et al., 2014). Finally, sea cucumbers studied in Mexico presented a much longer life than Holothuria (H.) grisea, as demonstrated by the reconstructed relative ages (^{Poot-Salazar et al., 2014}). This difference is probably associated to specific characteristics of the genus Isostichopus, since Holothuria arguinensis in Portugal partially overlapped the relative ages (^{Olaya-Restrepo et al., 2017}) of the hereby study.

The sea cucumber Holothuria (H.) grisea is a slow moving and unprotected invertebrate easily found in intertidal and subtidal zones along the Brazilian coast, therefore being susceptible to effortless manual collection during low tides or by snorkeling. It seems that in Santa Catarina these wild populations are still abundant but the threats of unregulated and unreported fishery are present. Therefore, urgent actions are required to protect these greatly important elements of marine ecosystems before the stocks become over-harvested or depleted, as reported for several species in various countries (^{Conand, 2017}; Conand, 2018; ^{Purcell et al., 2013}). Once the wild stocks have collapsed there is little capacity for natural recovery (^{Purcell et al., 2023}). We advocate that the development of sea cucumber aquaculture in Brazil is the only sustainable way to supply these highly prized organisms to the avid Asian markets while preventing wild stocks from being depleted.

Ethical statement: the authors declare that they all agree with this publication and made significant contributions; that there is no conflict of interest of any kind; and that we followed all pertinent ethical and legal procedures and requirements. All financial sources are fully and clearly stated in the acknowledgments section. A signed document has been filed in the journal archives.