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

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.

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

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:

Estimator of the total population size expressed in kg (Nȳ_st):

Estimator of the variance of the total size of the population $\widehat{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²_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:

with variance

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²_t is the variance of the total catch in the study area, m_t is the number of fishing operations and $\widehat{V}$(BT) 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 at (expressed in km) for each fishing operation resulted from

ai 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 $\overline{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²), 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 (${\overline{CPUA}}_i$ ). The total abundance of each i-th interval (N_i) was obtained by multiplying the (${\overline{CPUA}}_i$ ) 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, $\overline{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 $\widehat{SD}$ (BT) is the standard deviation of the total biomass. The abundance is calculated using Equation 10.

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₁= ${\overline{x}}^{-2}$ /S²–$\overline{x}$, subject to the following conditions: K/$\overline{x}$ >6 if $\overline{x}$ is low, K>13 if $\overline{x}$ is high or

$\left(\frac{k(k+\overline{x})(k+2)}{\overline{x}}\right)\ge 15$ if $\overline{x}$ is moderate.

If none of these conditions can be verified, then it is considered that K₁ is not suitable, and the parameter is estimated by means of the following equation:

For any case, p=$\overline{x}$/K . (18)

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

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.

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.

The niche breadth per sampling station ranged between 60 and 940 km², with an average of 242 km², 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.

The value of the parameters k and p (k²=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).

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