Potential biomass and distribution of octopus in the eastern part of the Campeche Bank (Yucatán, Mexico) ; Biomasa potencial y distribución de pulpo en el este del banco de Campeche (Yucatán, México)

1 Centro de Investigación y de Estudios Avanzados del IPN, CP 97310, Mérida, Yucatán, Mexico. (OA) E-mail: otilio.avendano@unicach.mx. ORCID iD: https://orcid.org/0000-0001-7644-8993 (IV-A) (Corresponding author) E-mail: jvelazquez@cinvestav.mx. ORCID iD: https://orcid.org/0000-0003-3216-2007 2 Universidad de Ciencias y Artes de Chiapas, Tonala, Chiapas, C.P. 30500, Mexico. 3Universidad Marista de Mérida, Mérida, Yucatán, C.P. 97300, Mexico. (AH-F) E-mail: ahernandez@marista.edu.mx. ORCID iD: https://orcid.org/0000-0003-1900-9868 (AC-J) E-mail: acuevas@marista.edu.mx. ORCID iD: https://orcid.org/0000-0002-8230-5021 4 Facultad de Ciencias Económicas y Empresariales, Universidad de Vigo, C.P. 36310, Vigo, Spain. (CF-J) E-mail: cjardon@uvigo.es. ORCID iD: https://orcid.org/0000-0003-0888-2055 5ECOBIOMAR, Instituto de Investigaciones Marinas, CSIC, C.P. 36208, Vigo, Spain. (AG) E-mail: angelguerra@iim.csic.es. ORCID iD: https://orcid.org/0000-0001-6716-3646


INTRODUCTION
The family Octopodidae contains more than 200 valid species. However, global octopus fisheries are based on only around 30 species, which are harvested in a range of fisheries from subsistence catches through to valuable, large-scale commercial fisheries (Norman et al. 2014). Since some of these catches are never properly identified, they are commonly reported as Octopus vulgaris. Due to the new findings in octopus taxonomy (Norman et al. 2014, Amor et al. 2017, Van Nieuwenhove et al. 2019, all landings reported under that name should be reviewed (e.g. Solís-Ramírez 1994, Norman et al. 2014, Emery et al. 2016). This confusion occurs in Mexico, one of the world's main producers of octopus, where two species are reported to be caught, O. maya and O. "vulgaris", the latter accounting for between 30% and 50% of the total landings ( Fig. 1) (Velázquez-Abunader et al. 2013, Sauer et al. 2020. For many years, managers and researchers have assumed that O. maya and O. "vulgaris" occupy different areas on the continental shelf of Yucatán, separated by depth (Salas et al. 2008, Beléndez-Moreno et al. 2014, DOF 2016. For instance, areas of higher abundance of O. maya were described along the coast, decreasing towards deeper waters (30 m deep at approx. 88°W) and increasing towards the west of Campeche Bank (Solís-Ramírez and Chávez 1986, Gamboa-Álvarez et al. 2015, Avendaño et al. 2019. Despite its importance, there has been no stock assessment of O. "vulgaris" in the Yucatán Peninsula (Salas et al. 2008, Sauer et al. 2020, which is the largest octopus-producing region in Mexico. In addition, some authors have noted that O. "vulgaris" from Yucatán has been misidentified (Amor et al. 2017, Ritschard et al. 2019, Van Nieuwenhove et al. 2019). It will therefore be referred to as O. "vulgaris" Type II in this article.
The National Institute of Fisheries of Mexico (IN-APESCA) has focused only on O. maya to determine management measures due to the difficulties of monitoring the fishing zones away from the coast where O. "vulgaris" Type II occurs (>30 m depth) (Gamboa-Álvarez et al. 2015, Ávila-Poveda et al. 2016. The assessment of O. maya included a total allowable catch (Salas et al. 2008, Velázquez-Abunader et al. 2013, DOF 2016. Other regulations for that fishery include a minimum landing size, a closed season for fishing and authorized fishing gear (DOF 2016), none of which consider the characteristics of O. "vulgaris" Type II.
Octopus production in Yucatán amounts to about 20000 t per year ( Fig. 1) (Velázquez-Abunader et al. 2013), and it is considered that the O. maya catch has reached the maximum sustainable yield (Jurado-Molina 2010, Beléndez-Moreno et al. 2014. However, the production of O. "vulgaris" Type II has shown a tendency to increase over time (Fig. 1), and in recent years has accounted for almost 50% of the total catch (Velázquez-Abunader et al. 2013). The catch of O. "vulgaris" Type II has increased due to the expansion of fishing activities and the introduction of a new medium-scale fleet (vessels between 12 and 15 m in length), which operates in deeper zones (~30 m deep), where it overlaps with the large-scale fleet (>15 m in length) (Salas et al. 2019, Beléndez-Moreno et al. 2014. This new medium-scale fleet also overlaps with the artisanal fleet (boats of 8 to 12 m in length) that targets O. maya in shallow waters (DOF 2016, Salas et al. 2019. Conducting fisheries management with limited knowledge of the biology of a species is risky (Pauly 1984, Leporati et al. 2015, and in some cases may result in overexploitation (Norman et al. 2014, Lima et al. 2017. In this regard, the Fishery Management Plan for the Octopus of the Gulf of Mexico and the Caribbean Sea (DOF 2016) underlines the need to enhance our scientific knowledge of O. "vulgaris" Type II because of the risk posed by the continuous growth in its landings.
There are several techniques for stock assessment, and their use depends on the information available and species characteristics (Hilborn and Walters 1992, Nevárez-Martínez et al. 2000, Hernández-Flores et al. 2015. Given that existing long time series and most of the information used for stock assessment of octopus from Yucatán are based on O. maya, the objective of this study was to assess the biomass status of O. "vulgaris" Type II and to determine its distribution on the northeast continental shelf of Yucatán at depths between 30 and 60 m. The data for performing this assessment were obtained from surveys conducted during scientific research cruises carried out independently from the fishery, so this information provides a preliminary analysis of the dynamics of a resource that has been catalogued with a status of potential exploitation. The analysis also explores the area of overlap between O. "vulgaris" Type II and O. maya.

Study area
The study area is located in the northeastern zone of the Yucatán Peninsula, on the continental shelf known as Campeche Bank, between 30 and 60 m depth (Fig. 2). The zone is influenced by the Loop current, which enters the system through the Yucatán channel, circulating towards the Straits of Florida and generating seasonal upwellings from May to August that stratify the water column, forming colder zones at greater depth (Tester et al. 1991, Lee and Williams 1999, Paris et al. 2005. From October to February, it is influenced by cold and dry air masses coming from the north of the continent (Canada and United States of America), which, when merging with air masses of tropical origin, mainly from October to March, produce strong frontogenesis (gusts >100 km h) and temperature drops (Zavala-Hidalgo et al. 2002, Morey et al. 2006, Dubranna et al. 2011, locally known as nortes (Enríquez et al. 2010).

Field work
Four research cruises were conducted on a semi-industrial vessel, from the port of Progreso, Yucatán, before and after the fishing season (from 1 August to 15 December). Fishing operations targeted O. "vulgaris" Type II at 29 stations on average (±2) (Fig. 2). The cruises took place in five months: May-June and July 2016, and December 2016 and January 2017.
The method used to collect octopus samples was the traditional fishing technique known as "gareteo". This technique consists of placing one bamboo pole (jimba) in the bow and another in the stern of a small wooden boat (<4 m long), known locally as alijo. These bamboo poles were approximately 4 m long (for the medium-scale fleet) (Salas et al. 2008, Sauer et al. 2020, and each had two monofilament lines with a specimen of Diplectrum sp. or Haemulon sp. tied to the tip as bait. Five alijos were deployed per cruise and each one had a GPS to trace the route and thus determine the swept area. The effective fishing time was standardized in one sampled by three hours per alijo (five) at each sampling station; thus, initial and final time of the fishing operations were recorded. At each fishing station or site, the total number of organisms caught (N t ), total catch weight (W t ) and total weight of each organism (w i ) were recorded.

Determination of the area per sampling station
Thiessen polygons (Brassel and Reif 1979) were deployed to identify the spatial arrangement of the sample as a function of distance between the sampling points. This allowed us to determine the optimal (constant) area of representativeness of each sampling point to apply any of the three area-based methods: stratified random method, swept area and weighted swept area. This procedure allows us to assign a real significance to sampling point values and therefore tells us the area on which the octopus abundance calculation was made according to the captures at each sampling station (Avendaño et al. 2019). The Thiessen polygons and the area of each polygon were estimated using the ArcMap 9.2 software (Sawatzky et al. 2009).

Biomass estimates
Four statistical methods were applied: the stratified random method (Cochran 1978), the swept areas method of Rosenberg et al. (1990), the Geostatistics biomass model (Rivoirard et al. 2008) and the method of weighted swept areas (Avendaño et al. 2019). These methods increase the precision of biomass estimation, thus generating robust data concerning the availability of octopus in the deepest areas (30 to 60 m) in the eastern part of the platform of the Yucatán Peninsula.
The stratified random method (Cochran 1978) uses the frequency of the distribution of the total weight of the catch per station, which must be classified by strata. The formula for calculating the biomass by means of the stratified method was obtained from the average count (kg) in the i-th stratum i(ȳ i ): Estimator of the variance of ȳ i : (2) Estimator of the total population size expressed in kg (Nȳ st ): ( Estimator of the variance of the total size of the population V(Nȳ st ): For the estimation of 95% confidence intervals: where N i is the total number of units sampled in the i-th stratum, L is the number of strata, n i is the number of units sampled in the sample of the i-th stratum, ȳ i is the average weight in the i-th stratum, and S 2 i is the variance of the counts in the i-th stratum.
The method of swept areas of Rosenberg et al. (1990) considered the catch in weight (biomass) obtained from the swept area, dividing the capture zone according to the observed abundances of the organisms: where Y t is the total catch in the study area, A t is the total area of the study zone, a t is the total area swept by the five alijos, S 2 t is the variance of the total catch in the study area, m t is the number of fishing operations and V(B T ) is the variance of the total biomass. In this case, a i represented the swept area covered by each alijo, so the total swept area a t (expressed in km) for each fishing operation resulted from a i was calculated using the following equation: where D i is the distance travelled by the i-th alijo, estimated using the tracking data recorded by the GPS, and LJ i is the length of the bamboo rods or jimbas installed in the i-th alijo (LJ i =8 m). Finally, the total abundance was calculated using the following equation: where TW is the average weight of the octopus obtained in biological samples. The biomass calculated by the geostatistical method was determined by previous calculation of the catch per unit area (CPUA, expressed in number of octopus per km 2 ), obtained by dividing the number of octopus captured by the area corresponding to each sampling station. The spatial correlation of the CPUA was determined by means of omnidirectional empirical variograms (isotropic), which measure the correlation between the semi-variance generated by all the differences between the data pairs separated by a previously established equal interval of distance (h) (Hernández-Flores et al. 2015). Later, a Kriging interpolation analysis was applied, giving the densities through interpolation nodes by the method of nearest neighbours (Cressie 1992), and generating a spatial structure that depends on the spatial arrangement of the population (Webster and Oliver 2007).
The empirical variograms were obtained using the following equation: where γ(h) is the value of the variance for each distance h, N(h) is the number of pairs of observations separated by h, C(x i ) is the CPUA observed at site x i and C(x i +h) is the CPUA observed at any other site separated by a distance h for each site x i . The obtained interpolations were divided into CPUA intervals, obtaining an average value for the i-th interval (CPUAi). The total abundance of each i-th interval (N i ) was obtained by multiplying the (CPUAi) by the total area covered per i-th interval; therefore, the total abundance (N T ) was obtained using the equation and the biomass was obtained using the equation The method of weighted swept areas proposed by Avendaño et al. (2019) analyses the density of organisms as a specific datum of the area represented by each station. That is, the abundance at a station is representative of the area of that station. Later, the total biomass is calculated from the sum of the biomasses obtained in each sampled area using the following formula: with a standard deviation of where Y i is the total of the catch in the i-th stratum, Y is the average of the catches, A i is the total area in the i-th stratum, a i is the area swept in the i-th stratum and  SD (B T ) is the standard deviation of the total biomass. The abundance is calculated using Equation 10.

Spatial distribution
The equation modified by Avendaño et al. (2019) from the method described by Guerra (1981) was used to determine the type of distribution, calculating the average probability of occurrence of octopus per fishing station, as well as, the type of distribution, and estimating the parameters of the negative binomial distribution p and k: To determine whether the type of distribution of the octopus is homogeneous or aggregated, a random distribution was developed, assuming a negative binomial distribution (k), calculated as K 1 =x -2 /S 2 -x -, subject to the following conditions: If none of these conditions can be verified, then it is considered that K 1 is not suitable, and the parameter is estimated by means of the following equation: For any case, p=x -/K .
Once the parameters had been estimated, a comparison was made between the distribution function of the total of samples and the theoretical distribution assuming a negative binomial type distribution by means of a goodness-of-fit test to check whether there was in fact a grouping (Zar 1999).

Spatial distribution of abundances
Empirical distribution models were generated using the variograms estimated in the geostatistical analysis because there was no previous evidence about the distribution of O. "vulgaris" Type II on the continental shelf of Yucatán. The model that best described the observed data was selected using two criteria: 1) the model that most closely explained the spatial distribution of the observed data, and 2) the minimum value of the sum of squares of the residuals, as well as the spatial distribution of the variances.

Analysis of geographic overlap between species
The distribution of O. "vulgaris" Type II was compared with the distribution of O. maya obtained by Avendaño et al. (2019) during the same period in the same study area. This analysis should corroborate the fact that their distribution zones are separated by depth, as described elsewhere (DOF 2016). The geographic overlap (D) for both species was determined using the Schoener equation (Schoener 1968, Cornell 2011, Afkhami et al. 2014). The Schoener index computes the percentage of similarity of two species and estimates the range of overlap in the areas where both species occur, based on the abundance of each one. Therefore, the over-geographical space is defined with a range between 0 (no overlap) and 1 (full overlap), where p xi or p yi denotes the proportion assigned to the species x (or y) in cell i.
The total geographic superposition of the area (D) was calculated using the equation: For both cases: where x i and y i are the number of organisms of the species x or y respectively at the sampling station i. X and Y are the total number of observations per species x or y.

Potential biomass analysis
The niche breadth per sampling station ranged between 60 and 940 km 2 , with an average of 242 km 2 , according to the Thiessen polygons. The lowest average biomass was estimated during the May-June cruise (37.8±3.36 t) while the highest was estimated in December (189.56±11.6 t) ( Table 1). The coefficient of variation (CV expressed as a percentage) obtained for the biomass estimated by the four methods was lowest during the months of May-June and July (CV=13.3% and 13.9%, respectively), corresponding to the period before the fishing season, while biomass increased in December and January (CV=20.6% and 20.1%) after closure of the season.
Among the four models, the lowest potential biomass was calculated using the geostatistical model, with values 18.5% to 36.7% lower than those obtained using the other methods. Similar patterns of variation were observed for density, but the differences from one month to another were not as evident as they were for biomass. Thus, the highest average density was recorded in December 2016 (46.5 octopuses per km -2 ), and the lowest average density was observed during the May-June 2016 cruise (6.8 octopuses per km -2 ). Finally, the weighted biomass method produced the second lowest values of biomass and abundance (Table  1), except during the May-June 2016 cruise.

Analysis of the distribution pattern
The value of the parameters k and p (k 2 =8, p=0.7) of the negative binomial distribution showed that individuals of O. "vulgaris" Type II tended to aggregate in patches. This species occupied the entire area, with density gradually increasing to a maximum in December and then decreasing to January. The variogram (Fig.  3) shows that as the distance between the sampling points increases, the O. "vulgaris" Type II CPUA (octopus km -2 ) increases towards shallow waters east of the Campeche Bank. This is due to the distribution of the species (Fig. 4).
Denser aggregations of O. "vulgaris" Type II were recorded towards the southeast of the Yucatán shelf. Cruises that took place before the fishing season (May-June and July 2016) showed the formation of high-density patches to the southeast of the study area (Fig. 4).
At the end of the fishing season (December 2016 and January 2017), the octopuses tended to be distributed heterogeneously over the sampled area and, although patches continued to form, they were less extensive (Fig. 4).

Geographical overlay analysis
In order to identify the overlapping geographic areas of O. maya and O. "vulgaris" Type II, a distribution analysis of both species was conducted using data collected during the same period (Fig. 4). The distribution of O. maya was obtained from Avendaño et al. (2019), and the distribution of O. "vulgaris" Type II was determined in the current study. The geographical overlap according to the Schoener index (1968) was relatively low for the May-June and July samplings (D=0.30 and D=0.39, respectively), which corresponded to the period prior to the fishing season. These values increased after the end of the fishing season, with higher values of overlap in December 2016 and January 2017 (D=0.60 and D=0.62, respectively).  An overlap in the distribution of the two species was observed during all the research cruises. This was more evident towards the south of the study area, where medium to low overlap was observed (D i ≈0.34) during the months prior to the opening of the fishing season; however, at those sites, O. maya was dominant in the west (approximately 88° W), while O. "vulgaris" Type II was dominant in the east (approximately 87° W). In the months following the fishing season, the heterogeneous distribution of smaller patches of both species caused the degree of overlap to increase throughout the study area (D i ≈0.61).

DISCUSSION
A total allowable catch, which is a common measure in cephalopod fisheries, is calculated in the octopus fishery of the Mexican Atlantic coast (Jurado-Molina 2010, Norman et al. 2014). Although this measure is not respected, it is aimed at O. maya (Salas et al. 2008, Beléndez-Moreno et al. 2014, DOF 2016, while O. "vulgaris" Type II can be freely caught. This study analyses for the first time the distribution and relative abundance of O. "vulgaris" Type II on the continental shelf of the Yucatán Peninsula between 30 and 60 m depth. This analysis was carried out through systematic sampling, independent from fishery activities, covering the potential fishing area. Because any precursor of fisheries management measures must be based primarily on the identification of resource distribution (Pierce and Guerra 1994), this analysis provides an instantaneous estimation of relative abundance in that area in four different periods (Rosenberg et al. 1990, Nevárez-Martínez et al. 2000.
The results on the biomass calculated using the four methods were relatively consistent (CV<21%), and discrepancies found between the methods could be due to the differences in the representativeness of the area of influence assumed by each one (Pierce and Guerra 1994, Nevárez-Martínez et al. 2000, Hernández-Flores et al. 2015. For example, in three of the methods (the swept area method, the weighted swept area method and the geostatistical biomass method), the distance between sampling stations is a key element in the calculation. In contrast, in the stratified method, the number of intervals is the key element, so for this method it is necessary to apply a stratification or classification from the start of the analysis (Nevárez-Martínez et al. 2000, Avendaño et al. 2019. The precautionary approach recommends choosing the lowest value calculated from the models to determine the octopus abundance in the area of study for any potential management measure, which in this case was obtained using the geostatistical method. As expected, our results showed that changes in the biomass and distribution of O. "vulgaris" Type II oc- curred before and after the fishing season in 2016. The increase in the relative biomass of O. "vulgaris" Type II from June to December (47%) and the change in its distribution pattern in those periods (Fig. 4) indicated that the stock reached its maximum abundance at the end of the fishing season (Table 1). This pattern differed from that observed for O. maya, which reached its highest value in October and gradually decreased to its lowest value in December (Table 2). This could be related to the fishing effort dynamics, as most fishermen target O. maya at the start of the season, but as its abundance decreases, they shift to O. "vulgaris" Type II. The commercial landings of O. "vulgaris" during the 2016 season peaked in December (1865 t), which coincided with the highest relative biomass estimated in this study using data collected independently from the fishery. To explain the catch levels of O. "vulgaris" recorded by the commercial fleet, the instantaneous abundance calculated in this study (average 189.6 t) for December 2016 could be 12 or 15 times smaller than the value of the relative abundance in the entire fishing area, meaning that the biomass in that month would be in the order of magnitude of between 2300 and 2800 t.
The propensity for O. "vulgaris" Type II to form aggregations, as observed in this study, infers that the behaviour of the species in this area is similar to that observed in O. "vulgaris" Type III in South Africa (Oosthuizen and Smale 2003) and O. insularis in Brazil (Leite et al. 2009). The higher abundance of O. "vulgaris" Type II towards the southeast zone and the formation of zones of aggregation in the months of May-June and July, as well as the heterogeneous distribution for the months of December and January, could be related to changes in environmental conditions, as suggested for other coastal cephalopods (Pecl and Jackson 2008, Pierce et al. 2008, Gamboa-Álvarez et al. 2015. This hypothesis stems from evidence that the Loop current and the Yucatán current have a strong influence on some other resources in the area (i.e. Epinephelus morio and Octopus maya) (López-Rocha et al. 2013, Gamboa-Álvarez et al. 2015. It is also possible that O. "vulgaris" Type II paralarvae are recruited in the eastern part of the Yucatán shelf due to the oceanographic conditions in that area (i.e. upwelling and nutrients [Enríquez et al. 2010]).
In the current study the northeast of the Yucatán Peninsula showed a geographical overlap of both species, ranging from low to medium overlap; this area was previously considered exclusively as an O. "vulgaris" Type II area. It is possible that this species is in fact more abundant in deeper waters, since it has not been recorded in the landings of artisanal fishermen. However, O. maya can occupy deeper waters than previously thought (Avendaño et al. 2019). Our results show that there are no clear limits on the distribution of the two species in waters at depths greater than 30 m, with both species being able to occupy the same space. The dominance of one or other of the species will depend on intrinsic biological characteristics and the prevailing ecological conditions, as well as on competition for space (e.g. success in recruitment or faster growth) (Leite et al. 2009, Pierce et al. 2008, Rodhouse et al. 2014). This scenario is consistent with the fishers' and researchers' idea that O. maya and O. "vulgaris" Type II occupy different habitats on the continental shelf, according to depth (Salas et al. 2008, DOF 2016. Thus, most abundant fishing grounds for O. maya were located along the coast, decreasing towards deeper waters (30 m depth at approximately 88°W) and increasing towards the west of the peninsula (Solís- Ramírez and Chávez 1986, Gamboa-Álvarez et al. 2015, Avendaño et al. 2019.
In conclusion, our results showed the changes in the distribution and abundance patterns experienced by O. maya and O. "vulgaris" Type II over a short period of six months. The results suggest that O. "vulgaris" Type II occupied the entire sampled area and formed patches, especially during December and January. Moreover, our findings revealed a geographic overlap of both species, whose amplitude changed throughout the year according to the geographic position: O. maya dominated at approximately 88°W, while O. vulgaris Type II dominated towards the southeast at 87°W.
Although in general the total abundance of both species was similar in the entire area of study, Octopus "vulgaris" Type II showed higher abundance at depths of less than 40 m. Consequently, it would be necessary to calculate the maximum sustainable yield of this species to develop appropriate management measures. The two species must be managed jointly as they overlap in great part of the fishing grounds. We also propose that it would be advisable to identify and create protected areas, which would provide refuge for both species at depths greater than 40 m. These areas would protect O. maya and O. "vulgaris" Type II adults, resulting in improved reproductive success, as mentioned by Velázquez-Abunader et al. (2013) and Avendaño et al. (2019). also thank Miguel Ángel Cabrera, Jesús Miguel Soto Vázquez and Luis Enrique Ángeles-González for their support in the field work.

Financial support
The results of this study were obtained as part of the project "Distribución, reproducción, biomasa y