Overexploitation and More Than a Decade of Failed Management Leads to No Recovery of the Galápagos Sea Cucumber Fishery

Isostichopus fuscus is the most important sea cucumber species exploited in the Eastern Tropical Pacific. It was the most important fishery in the Galápagos in the early 2000s until overfishing led to its collapse and a 5-year total fishing ban was established (2016–2021) to try to recover it. The management of I. fuscus in the Galápagos has always considered the population density as an indicator for decision-making. The objective of this study is to review the density as an indicator of the population status of I. fuscus using a stock-production model incorporating covariates methodology. For the first time, population and fishing parameters (K, r, and q), reference points (MSY, BMSY, FMSY, and DMSY) and indicators (B/BMSY and F/FMSY) were estimated for I. fuscus. The results indicate that the management measures have not prevented the overexploitation of this species for more than a decade. The goal of the I. fuscus management plan in the Galápagos, i.e., the recovering the fishery in a non-fishing scenario, will not be met by 2030. To accomplish its recovery six recommendations are proposed, including to extend the total ban of the fishery and to change the current management indicators to B/BMSY and F/FMSY. This study evidences that management measures taken with little scientific basis can have a pervasive effect on natural resources.


INTRODUCTION
Sea cucumbers have a high demand and value in markets of Asia because they are considered a food delicacy and have an important cultural in traditional medicine (Purcell, 2010;Fabinyi, 2012;To and Shea, 2012;Purcell et al., 2018). At the same time, sea cucumbers are very vulnerable to overexploitation due to the ease of capture that contributes to accelerated exploitation, increased demand, inefficient fisheries management, low recruitment, high longevity and density-dependence reproduction (Kinch et al., 2008;Purcell, 2010). Several countries and territories have established moratoria on sea cucumber fishing due to population declines (e.g., Australia, Mauritius, Mayotte, Papua New Guinea, Salomon Islands, Venezuela, and Ecuador including Galápagos) (Purcell et al., 2012;Dirección del Parque Nacional Galápagos et al., 2016). The International Union for Conservation of Nature has classified nine species of sea cucumbers as vulnerable and seven as endangered, the latter includes Isostichopus fuscus (IUCN, 2020). The Convention on International Trade in Endangered Species has three species of sea cucumbers listed on Appendix II and the only species listed on Appendix III is I. fuscus (CITES, 2020).
It is common for density-based reference points to be used to determine the population status of sea cucumbers. Purcell et al. (2009) considered reference points for commercial sea cucumbers in New Caledonia in a tentative and subjective way. These authors considered densities <1 sea cucumber/100 m 2 as low and densities < 0.3 sea cucumber/100 m 2 to be at, or near, a critical level at which populations will fail to repopulate effectively. The authors suggested to interpret these reference points in a meta-population dynamics context and that the suggested thresholds will be more realistic for some species than others.
The Secretariat Pacific Community developed reference densities for 18 species of sea cucumbers species of the centralwestern Pacific Ocean (Pakoa et al., 2014a). These densities are an average of the 25% of the highest densities recorded for each species from 2002 to 2012 in 91 sites of 17 countries. These reference densities were proposed as baseline to check for healthy abundances in a "rule of thumb" way in absence of site-specific density reference for sea cucumbers (Pakoa et al., 2014a). Pinca et al. (2010) proposed to use densities observed to infer at which one a stock is considered depleted or in healthy condition. These authors suggested reference points for Holothuria whitmaei, Bohadschia argus, and Holothuria atra using their densities collected in Pacific Island countries and territories. In Seychelles, Aumeeruddy et al. (2005) used densities to determine the status (i.e., underexploited, unexploited, exploited, and overexploited) of more than 20 species of sea cucumbers.
In few monospecific sea cucumber fisheries, density is used as an indicator to open or close fisheries. Examples are Cuba for Isostichopus badionotus (Hernandez-Betancourt et al., 2018), Canada for Parastichopus californicus (Fisheries and Oceans Canada, 2019) and Galápagos, Ecuador for I. fuscus.
Surplus production models have been applied to know the population status of three species of sea cucumbers previously, in addition to I. fuscus in the present study. Hernandez-Betancourt et al. (2018) applied the dynamic surplus production model of Schaefer to assess the population of I. badionotus in the southern coast of Camagüey, Cuba. Koike (2017) used the Pella-Tomlinson surplus production model to assess the population of Holothuria fuscogilva in Seychelles. Hajas et al. (2011) developed a surplus production model and compared it with the Schaefer model for Parastichopus californicus in British Columbia, Canada. Bradbury et al. (1996, 1998) used Schaefer model and Woodby et al. (1993 applied Caddy surplus production model for P. californicus in Washington and Alaska, United States, respectively. I. fuscus is the main commercially exploited sea cucumber in four countries of the Eastern Tropical Pacific: Mexico, Panama, Peru and Ecuador (Toral-Granda, 2008). Exploitation of I. fuscus in mainland Ecuador began in 1988 to supply Asian markets but soon decayed due to a collapse of the population from overfishing. Since 1992, there is a fishing ban in mainland Ecuador that prompted a migration of the activity to the Galápagos Islands (Carranza and Andrade, 1996).
In the Galápagos Marine Reserve (GMR), the I. fuscus fishery was the most important at the beginning of this century until it was considered overexploited in the mid-2000s (Dirección del Parque Nacional Galápagos, 2015. In 2016, a total closure for 5 years was established for this fishery in order to achieve its recovery in the GMR (Dirección del Parque Nacional Galápagos et al., 2016).
The management of the I. fuscus fishery in the GMR has always considered the density of its population as an indicator on which to decide to open or close the fishery since it has been assumed as a representation of the population status. In 2002, the GMR co-management system approved two reference points on which to decide the opening of the I. fuscus fishery each year: (1) 40 sea cucumbers/100 m 2 ; (2) That the catch per unit effort (CPUE) does not decrease for three consecutive years in three quarters of the macro-zones into which the I. fuscus fishing area in the GMR was divided 1 .
In 2008, these indicators were changed into a traffic light system in which the density in the west of Isabela Island indicates the population status of I. fuscus in the GMR. A density greater than or equal to 21 sea cucumbers/100 m 2 indicates a healthy population; if the density is equal to 11 or less than 20 sea cucumbers/100 m 2 the population is in a recovery status and; if it is less than 11 sea cucumbers/100 m 2 it is in a critical status. This system contains decision rules where, if the I. fuscus population is healthy or recovering, the fishery opens with 60 days of fishing, a total allowed catch (TAC) and closure of critical areas (e.g., recruitment areas or areas with low densities). If the population has a critical status, the fishery remains closed . This system is currently in place.
Given this historical framework, the objective of this study is to use a surplus production model to review the density as an indicator of the population status of I. fuscus and its use as a decision rule for the management of the fishery of this species in the Galápagos. This study is the first to estimate population and fishery parameters of I. fuscus, which will contribute to the knowledge and management of this species in the Tropical Eastern Pacific region.
This study is the second stage of a reengineering of the management of this resource in the Galápagos. The first stage was to estimate the age and growth, to update the total mortality and the fishing and habitat area of this species in the GMR in order to propose a change in the methodology to estimate the TAC. In fact, Ramírez-González et al. (2020) recommended the use of a dynamic surplus production model to estimate several reference points for I. fuscus in the GMR. The third stage will be to carry out spatially explicit analyses of the population of I. fuscus throughout the GMR and to determine the implications for its fishery management. This reengineering has the ultimate goal to bring scientific inputs to the GMR co-management system, before 2021, when negotiations to open or keep the sea cucumber fishery closed will be back on the table and decisions on the management of this fishery will be taken 2 , especially under the new scenario of the recent economic crisis due the COVID-19 pandemic.

Study Area
The GMR is located 1,240 km to the west of Guayaquil city, Ecuador and has its geographical center at 96 • 46'W and 0 • 05'S. The reserve covers the entire marine area within 40 nautical miles measured from the baselines of the Galápagos Archipelago and its inland waters. It comprises a total area of ∼138,000 km 2 (Figure 1; Dirección del Parque Nacional Galápagos, 2014).
The GMR is a marine protected area of multiple permitted uses, including small-scale fishing, tourism, research and conservation. These uses were spatially restricted to areas defined by a coastal zoning from 2000 to 2015, where the no take zones covered 1% of the GMR total area and 77% of its coast (Heylings et al., 2002;Moity, 2018a). In 2016, new zoning was established to include also open waters of the GMR, where the no take zones cover 32% of the GMR's total area (Ministerio del Ambiente, 2016). Edgar et al. (2004) divided the GMR in the following five bioregions: Far northern, Northern, Western, Elizabeth and Central-southeastern (Figure 1). Each bioregion was defined on the basis of multivariate analysis of shallow subtidal rocky reef fish and macroinvertebrate communities. Each bioregion is characterized by mix of species derived from Indo-Pacific, Panamic, Peruvian, and endemic source areas. Both monitoring programs were designed to estimate the exploitable abundance of I. fuscus and were carried out in the fishing and habitat area of I. fuscus in the GMR (Figure 1). To estimate density, all the sea cucumbers inside a circular plot of 100 m 2 (radius of 5.64 m) were counted. In total, 6,887 plots were settled between <1 and 38 m depth (11.7 m mean; 6.1 m SD) (Supplementary Table 1).

Data Sources
Fishing data of I. fuscus from the GNPD datasets for years 1999 to 2005, 2011, and 2015 were also used. Datasets for years 2007 and 2008 were excluded because they did not contain enough variables for the CPUE standardization (see section "I. fuscus Standardized CPUE"). For the rest of the years there were no 2 Dirección del Parque Nacional Galápagos. Meeting minute, June 5, 2020. fishing data due to the sea cucumber fishing ban. For this study, only fishing data from landing ports monitoring were used because sampling with captain's logbooks and onboard observers where conducted only in 2000-2001 and 1999-2005, respectively. The landing ports monitoring consisted of persons who recorded the date and counted all the sea cucumbers from each arriving fishing vessel. The captain of each fishing vessel was also interviewed about the number of fishing days, vessel type, whether the vessel was towed to a mother boat, the number of divers, the fishing site, the type of diving and the fishing depth (in meters). For this study, only hookah diving as the type of diving was used because scuba diving and free diving had only 41 observations in total (vs. 11,438 observations for hookah diving).

Indicators and Reference Points
A stock-production model incorporating covariates (ASPIC) methodology (Prager, 1994(Prager, , 2016 was employed to estimate indicators and reference points for I. fuscus fishery in the Galápagos. ASPIC fits in a non-equilibrium way the following surplus-production models: the logistic production model of Schaefer (1954Schaefer ( , 1957 and Pella (1967); the generalized model of Pella and Tomlinson (1969) as reparametrized by Fletcher (1978) and; the Fox exponential yield model (Fox, 1970). The fit in a non-equilibrium way is important for sea cucumbers because it has been documented that several species do not have an annual recruitment, which violates the equilibrium assumption (Purcell, 2010).
ASPIC's outputs includes the following estimations with confidence intervals: carrying capacity (K), catchability coefficient (q), maximum sustainable yield (MSY), fishing mortality rate at MSY (F MSY ), and abundance at MSY (B MSY ). ASPIC also estimates projections of population abundance and fishing mortality rates.

I. fuscus Abundance
The annual abundances were estimated with the following equation modified from Wolff et al. (2012b) for I. fuscus in the GMR: Where B t is the abundance in number of sea cucumbers for the year t, d tj is the median of the density given in sea cucumbers/100 m 2 for the year t and monitoring zone j (obtained from the population monitoring programs) and A j is the fishing and habitat area of I. fuscus in the monitoring zone j. Because, in every case, the densities did not have a normal distribution (Shapiro-Wilk; p < 0.001), and were skewed toward zero, the median was used as a representative measure of central tendency of d tj . In 2012, no density sampling was done on the islands of San Cristóbal and Española, and in 2014, no density sampling was done in the Elizabeth bioregion; so the density values for these years and monitoring zones were assumed by averaging the previous year's density value with the following year's density value (Supplementary Table 2).
Since it is known that sea cucumbers can have spatial density patterns that depend on different environmental characteristics (Domíngez-Godino and González-Wangüemert, 2020), a Kruskal-Wallis test with 95% confidence was performed to see differences between the densities across the different monitoring zones of I. fuscus. Each monitoring site was assigned to a bioregion on the basis of spatial analysis using Moity (2019) bioregion dataset in ArcGIS 10.6.1. There were statistical differences between the three bioregions in which the monitoring programs were carried out (i.e., Western, Elizabeth and Central-southeastern). Because the Central-southeastern bioregion included more sampling zones in a wider spatial range than the other two bioregions, another Kruskal-Wallis test was conducted to see differences between the islands within this bioregion. There were statistical differences between Española Island and the rest of islands and between San Cristóbal island and the rest of the islands (Supplementary Table 3). With these results, the following zones (A j ) were determined for the annual abundance estimate: Western, Elizabeth, Central-SE Santa Cruz-Floreana, Central-SE San Cristóbal and Central-SE Española (Figure 1). The fishing and habitat areas of I. fuscus of each monitoring zone were calculated using the methodology of Ramírez-González et al. (2020) and the dataset of Moity and Ramírez-González (2019) (Supplementary Table 3).

I. fuscus Standardized CPUE
The CPUE was defined as the number of sea cucumbers caught per diver per day. Zero values of CPUE (n = 6) and CPUE in the Northern bioregion (n = 1) were not used for the standardization. A Linear Multiple Regression model was used to standardize the CPUE for each year. Prior the CPUE standardization, Generalized Variance-Inflated Factor (GVIF) tests were conducted to select the covariates that were not correlated. The covariates tested were year, landing port (Puerto Ayora, Puerto Baquerizo Moreno, Puerto Villamil), month, towed (yes, no), vessel type (fiberglass, wood), bioregion of fishing site (according to Edgar et al., 2004; and fishing depth. According to Zuur et al. (2010) covariates with GVIF values >3 are correlated. We sequentially dropped collinear covariates following Zuur et al. (2010)'s criteria until all covariates had GVIF values <3. The final explanatory covariates chosen were year as factor, towed as factor, bioregion as factor (Western and Elizabeth bioregions were pooled) and fishing depth as continuous (Supplementary Table 4). The tests were conducted The CPUE data was log-normalized. Normality of log-CPUE was tested visually by means of Q-Q plots and accumulative density plot (Supplementary Figure 1). Q-Q plots revealed deviance from normality, however, we kept using log-CPUE in the model since there are known issues with null-hypothesis significance test for non-normality with very large sample sizes, as it is the case in this study (n = 11,431). See Supplementary Material for the model validation.

ASPIC Fits
The ASPIC 7 software (Prager, 1994(Prager, , 2016 in BOT program mode was used to fit de time series with the logistic Schaefer model and to compute bootstrapped confidence intervals. The estimation method chosen was maximum a posteriori (i.e., a form of penalized likelihood using Bayesian priors) with 20,000 Monte Carlo searches. We performed 1,000 bootstraps with 95% confidence intervals as recommended by Prager (2016). The seed values of abundance at the start of the first year, K ratio (B 1 /K), MSY, F MSY and q were 0.5, 5.27 × 10 6 , 0.058 and 6 × 10 −6 , respectively. B 1 /K was determined as Prager (2016) suggests. MSY, F MSY , and q were determined with preliminary maximum likelihood fits of the model (See Supplementary Material for ASPIC 7 script). The contrast index and the nearness index were used to determine the goodness of fit of the model. The contrast index is the mean of B coverage proportions > and < B MSY , where a value = 0.5 is a good fit and value = 1 the best fit (Prager, 2011). The nearness index is the proportional closeness of any B to B MSY , it takes values from 0 to 1, where 1 is the best fit (Prager, 1994(Prager, , 2016. With the estimations of MSY, K, and B MSY , the intrinsic rate of population growth (r) was calculated with r = MSY/(K/4) and the

RESULTS
I. fuscus density range of all the plots analyzed was 0-860 sea cucumbers/100 m 2 from 2000 to 2019 (5 median; 2 25thpercentile; 9 75th-percentile). The CPUE range was 1-9,320 sea cucumbers/diver per fishing day (256 mean; 309 SD). Supplementary Table 1 shows the 2000-2019 time series of abundance estimated, CPUE standardized and total catches recorded of I. fuscus in the GMR.
The CPUE model fitted well to the CPUE standardized, meanwhile the abundance model fitted well as of 2003 (Figure 2). Schaefer's model had a good fit in relation to B and B MSY , the contrast index was 0.6 and the nearness index was 1. Table 1 shows the estimates of the population parameters and reference points of I. fuscus in the GMR.
According to the B/B MSY indicator, I. fuscus has been overexploited since 2004 and the value of F/F MSY has always been greater than 1. The 10-year projection of the B/B MSY indicator under a non-fishing scenario estimated that the I. fuscus population will not reach the B MSY in 2030 (Figure 3).

DISCUSSION
The use of density as reference point for the management of the I. fuscus fishery in the GMR has failed because it did not prevent the overexploitation of this species. The density estimates as an indicator have been based on observations without further verification and have not considered the spatial dynamics of this species. The density of 40 sea cucumbers/100 m 2 proposed as a "sustainability coefficient" of I. fuscus in the GMR in 2002, was based on observations where the reproductive aggregations of this species had at least this density 3 . Castrejón et al. (2007) said that this density was poorly estimated, since, until then, the average density of this species was never higher than 35 sea cucumbers/100 m 2 in the GMR. They noted also that it was not clear whether this reference point referred to the total fishing area or to a specific fishing area. Moreover, it has not been proven whether 40 cucumbers/100 m 2 is the optimal density for the reproduction of I. fuscus in the GMR.
The density of 11 cucumbers/100 m 2 proposed as indicator of critical status of I. fuscus population in the GMR, was based on the observation that the lowest catches of this species occurred with densities lower than this value on the west of Isabela Island . No rigorous analysis has been done to see the relation between catches and density of this species. Furthermore, as demonstrated in this study and in Glockner-Fagetti and Benítez-Villalobos (2016), I. fuscus has a differentiated spatial abundance, so the density of an area (e.g., west of Isabela Island) does not represent the density throughout its distribution range. Therefore, it is not recommendable to determine a single density value as an indicator of population status of I. fuscus in the GMR.
The D MSY estimated in this study (0.57 sea cucumbers/100 m 2 ) is quite close to the estimated for I. badoinotus in Cuba (0.45 sea cucumbers/100 m 2 ; Hernandez-Betancourt et al., 2018). However, it is much higher than the D MSY estimated for H. fuscogilva in Seychelles (0.002-0.005 sea cucumbers/100 m 2 ; Koike, 2017), although it is a much larger species (Purcell et al., 2012) that occurs naturally at lower densities (Pinca et al., 2010). The D MSY of I. badoinotus and H. fuscogilva are probably underestimated because the authors assumed that the entire study area is suitable habitat for each respective species. Our D MSY is probably also biased since the habitat area of I. fuscus was modeled with observations mostly in the south of the Northern bioregion and with a still not very accurate bathymetry for the GMR (Moity, 2018b).
Our D MSY is in line with the density-based reference points suggested by Purcell et al. (2009) for sea cucumbers in New Caledonia, in the sense that less than 0.3 sea cucumber/100 m 2 indicates a critical status. Of the 18 species with regional reference densities of the Secretariat of the Pacific Community, five exceed our D MSY : Holothuria atra (56 sea cucumbers/100 m 2 ), Stichopus chloronotus (35 sea cucumbers/100 m 2 ), Bohadischia similis (14 sea cucumbers/100 m 2 ) Holothuria coluber (11 sea cucumbers/100 m 2 ) and Holothuria scabra (0.7 sea cucumbers/100 m 2 ) (Pakoa et al., 2014a). It appears that I. fuscus requires relatively higher densities than many of tropical sea cucumber species to keep its population healthy.
I. fuscus has higher densities than those for other tropical sea cucumber species, considering the observed values. Pooling the annual densities of this study, the median density of I. fuscus per monitoring zone is nine, eight, five, four and six sea cucumbers/100 m 2 for Western, Elizabeth, Central-SE-Santa Cruz-Floreana, Central-SE-San Cristóbal, and Central-SE-Española, respectively. Of more than 25 species of tropical sea cucumbers studied in the Central Pacific (Purcell et al., 2009;Pinca et al., 2010;Pakoa et al., 2013Pakoa et al., , 2014b and Seychelles (Aumeeruddy et al., 2005;Koike, 2017), only H. scabra, Bohadschia vitensis, Stichopus spp., Holothuria hilla and H. atra had occasionally observed densities greater than four sea cucumbers/100 m 2 .
Observed densities of I. fuscus have varied with respect to its geographic distribution and time, with the highest densities recorded in its most northern (northwestern Mexico) and southern (Galápagos) distribution. In Baja California Peninsula densities of 30 sea cucumbers/100 m 2 of I. fuscus (Fajardo-León and Vélez, 1996), 27 sea cucumbers/100 m 2 , 15 sea cucumbers/100 m 2 (Glockner-Fagetti et al., 2016) and 2.8 sea cucumbers/100 m 2 (Reyes-Bonilla et al., 2008) have been recorded. In southwest coast of Mexico a density of 43 sea cucumbers/100 m 2 were recorded in 1991, years later less than one sea cucumber/100 m 2 (Nuño-Hermosillo, 2003) and 1.8 sea cucumbers/100 m 2 were recorded (Glockner-Fagetti and Benítez-Villalobos, 2016). In the El Pelado Marine Reserve (mainland Ecuador), a density range of less than one to three sea cucumbers/100 m 2 has been recorded for 2014-2015 (García, 2015). The maximum density ever observed of I. fuscus was 88 sea cucumbers/100 m 2 in the Western bioregion in the GMR in 2001.
It appears that I. fuscus is among the tropical sea cucumber species with the highest observed densities. Caution should be taken in interpreting these results as the density value may vary due to environmental factors, such as, ocean productivity (Reyes-Bonilla et al., 2008;, whether the population is or was subject to exploitation (Holguin-Quiñones et al., 2000), and time (e.g., month or day and night; . On the other hand, management measures related to controlling fishing mortality (e.g., TAC and fishing season) of I. fuscus in Galápagos did not prevent this mortality from being above the F MSY during the entire study period. From 2000 to 2008, the TAC was established using averages of catches from previous years. In 2011, the TAC was implemented by estimating MSY using the Cadima formula. In 2015, the TAC was issued at the discretion of the participants in the GMR co-management system since in this year the fishery should not have opened because it did not reach the reference point (Ramírez-González et al., 2020). The latter authors discussed the lack of rigor in the estimation of TAC in the first years and in 2015, they also concluded that the use of Cadima formula is no longer appropriate because its scope is very limited. All the TACs established for the I. fuscus fishery in the Galápagos were above our MSY estimate (Dirección del Parque Nacional Galápagos, 2015).
Illegal, underreported and unregulated (IUU) fishing of I. fuscus should also be considered. Schiller et al. (2013) found that the largest annual catch of this species occurred in 1994 (2,888 tons), when the fishery was not regulated in the Galápagos. They also calculated that 3,000 tons of this species were caught illegally between 1994 and 1999 in the Galápagos. In addition, in 2000, 2003, and 2015 the established TAC was exceeded (Reyes et al., 2013a,b;Ramírez-González et al., 2020). This IUU has most likely contributed to the overexploitation of I. fuscus in the Galápagos.
The K estimate of this study suggests that the I. fuscus population has suffered a strong impact in the GMR. Considering K as virgin abundance and with a conservative approach (using the K low confidence interval), the population of this species in 2019 was 33% of its virgin population. Wolff et al. (2012b) assumed that the abundance of I. fuscus in 2009 was already 50% below of its virgin abundance; our estimations indicate that in 2009 the abundance was 15% below (under a conservative approach), which is, in fact, the lowest percentage recorded during the study period.
The intrinsic growth rate estimated in this paper (r = 0.02) corresponds to that of a slow-growing population. We have only found two other r-values for sea cucumbers in the literature, in both cases our r-value is much lower than Cisneros-Mata (2016) for I. fuscus in Mexico (r = 0.5) and Hernandez-Betancourt et al. (2018) for I. badionotus in Cuba (r = 0.129). The difference with Cisneros-Mata (2016) is most likely, due to that they assumed the r-value from ecological and biological information on the species. The difference with Hernandez-Betancourt et al. (2018) could be due to the fact that the r-value of I. badionotus was estimated with a non-overexploited population and our r-value corresponds to an overexploited population. One thing that stands out is that, to our knowledge, this is the first time the population parameters of I. fuscus are estimated.
The impact of overexploitation led to the fact that, according to our 10-year projection, I. fuscus will not reach the MSY level in 2030 under a non-fishing scenario in the GMR. There are examples of sea cucumbers with slow recovery after a moratorium. Holothuria whitmaei still had depleted levels after a 7-year moratorium in Tonga (Friedman et al., 2011). H. scabra had a slow recovery after 6-year moratorium in Warrior Reef, Australia (Kinch et al., 2008).
However, there is an increasing trend in the I. fuscus population from 2016 onwards in the GMR. This implies that the impact of overexploitation did not exceed density thresholds where the reproductive potential of I. fuscus fails to recovery its population. The lowest median densities per monitoring zone were observed in 2009 (1-4 sea cucumbers/100 m 2 ). Therefore, probably this range is above the tipping point where I. fuscus can no longer recover in the GMR. This considering that the sea cucumbers are gonochoric (i.e., separate sexes) and have low mobility, which make them highly dependent on its density (Lovatelli et al., 2004).
The no recovery of I. fuscus in the short-term in the Galápagos is important because currently fishers and part of the local population have considered this resource as an alternative for economic reactivation in the face of the negative impacts of the COVID-19 pandemic. The results of this study contribute to realign the expectations about I. fuscus in the Galápagos and to start thinking that this resource is not going to be available in the long-term.
It is important to account for the climate variability as a factor that may have influenced our results. The higher abundances of the first years of the study (2000)(2001) may have been influenced by the 1998-1999 strong La Niña event, being one of the strongest after the 1988-1989 La Niña in the GMR (National Oceanic and Atmospheric Administration, 2013). The high productivity associated with such events (Palacios, 2004), especially on the west of Isabela Island and in Fernandina Island, where major fishing grounds occur, could have yielded a peak in I. fuscus reproduction and abundance, in the following years. There are studies that have demonstrated influence of climate variability on I. fuscus population in the west of the GMR (Hearn et al., 2005;Wolff et al., 2012a). This could result in a biased baseline in the model used in this study.
Absolut values of the parameters and reference points estimated should be taken with caution. Prager (2016) says that ASPIC estimate more precisely MSY, B/B MSY and F/F MSY and less precisely absolute levels of B, F, and q. He also mentions that B/B MSY and F/F MSY explain better the condition of the population, because in normalization, the estimate of q cancels out. Nevertheless, even considering the biases, there is evidence that the I. fuscus population in the Galápagos is overexploited and without a recovery to MSY levels in the short-term.
The biological objective for the I. fuscus fishery established in the Five-Year Fishing Calendar 2016-2021 (Fishing management plan of the GMR) states as "to recover the abundance of the sea cucumber populations ensuring a healthy population structure" (Dirección del Parque Nacional Galápagos et al., 2016). Our results show that the 5-year closure of this fishery will not meet this objective.
Finally, there are still research gaps to be addressed for the correct management of this sea cucumber species in the GMR. All the analyses have been done in the habitat and fishing area of I. fuscus, however, the species is widely distributed in the GMR (Ramírez-González and Reyes, 2018;Moity, 2018b) although it is not captured in the entire Reserve. The total abundance of I. fuscus in the GMR and the population spatial-temporal dynamics are unknown. It may be possible that the population of this species in non-fishing areas contribute to the recovery of the ones in the fishing areas due to a spill-over effect. In addition, we still do not know the potential impacts of climate change and climate variability on the I. fuscus population in the GMR. The impacts that the overexploitation of I. fuscus is causing on the ecosystems of the GMR are currently unknown. It is very important to conduct economic and social science research on this resource in the GMR and, to link it to biological and ecological research to understand the drivers resulting in its overexploitation.

Recommendations for the Management of I. fuscus Fishery in the GMR
The findings of this study suggest that the I. fuscus population is still not recovered from overfishing and the fishery should remain closed. Particularly, we recommend the following: 1. To use the trend in abundance of I. fuscus, considering the B MSY as limit reference point. The limit reference points are values that indicate an undesirable status of the fishery and should be avoided (Caddy and Mahon, 1995). 2. To open the I. fuscus fishery when its abundance in the GMR reach a positive trend above the B MSY value. 3. To restock I. fuscus in the GMR by collecting adult sea cucumbers present in the same area of a monitoring site and release them in the nearest no-take zone, according to GMR zoning. This intervention could be more costeffective than a captive-release program (Bell et al., 2008), but would need consideration of the legal, administrative and cost issues and to carry out control and surveillance in the no-take zones to ensure the success of this action. 4. To strengthen control and surveillance to reduce IUU fishing of sea cucumbers in the GMR to zero. 5. When the I. fuscus fishery opens in the GMR, to establish the TAC by multiplying the low confidence of interval of the F MSY by the B MSY . 6. To continue with the Sea Cucumber Population Monitoring Program of the GNPD and fishing sector of Galápagos. It is important to validate with empirical data the 10-years projection made in this study.

CONCLUSION
This is the first time that the population parameters of I. fuscus are estimated using a surplus production model. Under the categorization of Dowling et al. (2013), this fishery improved its data use and went from level 4 (empirical estimates of biomass based on fishery independent surveys) to level 2 (assessment of F and B based on fishery dependent and fishery independent data), where 0 is the maximum level. The management of I. fuscus have not prevented its overexploitation for more than a decade and have lead to absence of its recovery to the MSY level in the short-term in the GMR, even under a non-fishing scenario. This study evidences the impact that taking management measures with little scientific basis can have on natural resources.
To accomplish the objective of the GMR's fishing management plan of I. fuscus recovery, it is neccesary to continue the ban of its fishery and change its current management measures, such as the current indicator and reference points.

DATA AVAILABILITY STATEMENT
The datasets for this article are not publicly available because the Galápagos National Park is the owner and administrator of the data sets and decides on their availability upon request.
Requests to access the datasets should be directed to HR, hreyes@galapagos.gob.ec.

AUTHOR CONTRIBUTIONS
JR-G conceived the research idea, analyzed the data, and wrote the manuscript with inputs from NM and SA-V. HR was in charge of planning and conducting field trips, as well as collecting and storing field data. All authors contributed to the article and approved the submitted version.