Northwest Mexico is home to two of the country’s most productive fishery zones: the southern extension of the California Current Ecosystem, and the Gulf of California (GoC) (Cisneros-Mata, 2010, FAO, 2022). Indeed, from our own interrogation of the CONAPESCA database, we found that from 2006 to 2014, SSFs in Northwest Mexico contributed 52% of the total biomass generated by the nation’s marine SSFs (Ferrer et al., 2021).

Understanding the critical importance of SSFs in and around the GoC, colleagues with the GCMP and local small-scale fishers have worked collaboratively for over a decade to collect high-resolution fisheries monitoring data from in and around the Baja California Peninsula (see http://gulfprogram.ucsd.edu/slider-home/projects/). This monitoring program employs portable GPS tracking devices to populate a database with thousands of fishing tracks (> 25,000), about 5,000 of which are further appended with information related to catch and fuel consumption. These records within the GCMP database comprise the universe of data we used to generate the fuel intensity estimates described below.

A fuel intensity estimate (FIE) is a measure of fuel efficiency, telling us the fuel consumption per unit of X, where in this case X is one kilogram of wet weight ¹ , ² seafood. To generate FIEs, we divided the gasoline consumption (in liters) by the wet weight of catch reported for each single-species fishing trip identified/located in the GCMP database. To avoid the numerous uncertainties associated with partitioning fuel consumption among trips targeting multiple species (Vázquez-Rowe et al., 2012), we constrained our analysis to records for single-species trips. We also excluded any records that were obviously duplicated, incomplete, split, or erroneous (e.g., average boat speed > 80 km/hr). The resulting dataframe (n_records= 4,795) is included in the Supplementary Material section of this text, entitled “ Supporting Dataframe 1 ”. The raw data for this analysis (Mascareñas-Osorio et al., 2017) are available upon consultation at dataMares.ucsd.edu.

To generate stock assessment profiles (described in Section 2.4), we first needed to extract catch time series from the CONAPESCA database (Mascareñas-Osorio et al., 2018). The particular version of the CONAPESCA database that we used is a repository consisting of tens of thousands of records for small-scale and industrial fishing activities, targeting over 500 taxa, reported in all of Mexico’s Pacific states. Thus, to extract the relevant time series from this rather large database, we needed to identify the “Stocks of Interest” for which we would eventually require estimates of B/B_MSY (Section 2.4); we did so based on the definition of a “Stock,” and the two “Criteria for Inclusion,” we describe below.

We define a “Stock” (S) as any genus (G) living in fishing zone (Z), such that S = G_Z . Subsequently, we define a “Stock-Year” (S_Y ) as a stock (S) in a specific year (Y), such that one S_Y = G_Z,Y . So, for example, we consider Callinectes-Lower Pacific-2015 and Callinectes-Lower Pacific-2016 as two separate “Stock-Years” borne from the same “Stock” of swimming crab (Callinectes sp.). We defined fishing zones (Z) a priori, based on expert knowledge of the area and a large geographic separation of recorded ports ( Figure 2 ). The fishing zones relevant to our analyses are as follows: (1) the “Central Pacific” (CP), which includes ports circa Punta Eugenia; (2) the “Lower Pacific” (LP), which includes major ports on the Pacific side of Baja California Sur, near Bahía Magdalena; (3) the “Lower Gulf” (LG), which we’ve defined as the ports in and around La Paz, south of Loreto; and finally, (4) the “Upper Gulf” (UG), which spans the region from the Colorado River Delta to the Midriff Islands in the GoC. Note that, on average, small-scale fishers in the region only travel about 63 kilometers in a given fishing trip (see “ Supporting Dataframe 1 ”), and it is unlikely — though not impossible — that fishers in, for example, the LP would land and register their catch in the LG.

With these definitions in mind, we decided upon two “Criteria for Inclusion,” that is, two criteria used to determine whether or not a stock would be included in our final statistical analyses and therefore in need of assessment. First, we stipulated that a stock (S) can be found in the CONAPESCA database with ≥ 10-years’ worth of landings data; and second, that that stock has at least one Stock-Year (S_Y ) represented in the GCMP database with ≥ 10 paired observations of fuel consumption and catch. Ultimately, these “Criteria for Inclusion” produced a list of 19 “Stocks of Interest” for which we were able to extract landings time series from the CONAPESCA database, and which are associated with the 39 “Stock-Years of Interest” listed in Table S2 .

The reason for our first “Criteria of Inclusion,” in which we have stipulated that n years of CONAPESCA catch data must be ≥ 10, has to do with the number of years required to generate a reliable stock assessment. According to Froese et al. (2017), it is sometimes possible to generate an stock assessment profile with as little as five-year’s worth of catch data, however, the fewer the number of years, the larger the uncertainty. Moreover, for species with low or very-low intrinsic population growth (e.g., ocean whitefish Caulolatilus princeps), about 10 years’ worth of landings data are required to generate an informed prior about the fishery’s catchability and, in turn, its Maximum Sustainable Yield (Froese et al., 2017). The second criteria we have defined, that n paired observations of fuel use & catch in the GCMP database must be ≥ 10, has to do with the methods we used to test the relationship between B/B_MSY and fuel intensity among “Stock-Years of Interest” (described in Section 2.5); summarily, we wanted to include only those Stock-Years for which it would be possible to generate a reliable estimate of mean fuel intensity.

B/B_MSY is often used as a measure for stock status and a proxy for “overfishing,” where B is the current (available) stock biomass, and B_MSY is the estimated stock biomass at “Maximum Sustainable Yield” (MSY). MSY is the theoretical level of extraction equal to the rate of added population growth over an indefinite period of time — the highest level at which it is possible to extract fish from a population while still maintaining that stock’s standing biomass. By fishing a stock at MSY, the biological stability, economic value, and contribution to fishers’ livelihoods (derived from that stock) are more likely to be preserved in the long-term (Giron-Nava et al., 2021). Meanwhile, fishing a stock at levels above MSY is considered “overfishing,” which, if allowed to persist, will eventually diminish or crash that population (e.g., New England Cod; Pershing et al., 2015).

Herein lies a key distinction between “overfishing” and “overfished”. Overfishing, estimated as the current fishing mortality versus fishing mortality at MSY (F/F_MSY ), is the actual act of fishing beyond MSY at any given point in time. In contrast, B/B_MSY tells us if a stock is currently overfished, which, if it is, is usually the direct result of habitual or chronic overfishing. We contend that, in this respect, B/B_MSY carries some amount of “memory” vis-à-vis fisher behavior in response to stock status over time. We have therefore elected to use B/B_MSY as our proxy for overfishing, where biomass (B) exists on a gradient from “pristine” (B/B_MSY = 2) to “overfished” (B/B_MSY ≤ 0.5), such that when B is precisely equal to biomass at MSY, the ratio of B/B_MSY is equal to 1.

To generate B/B_MSY estimates for the 39 “Stock-Years of Interest” described above (in Section 2.3), we applied a data-poor method for fishery stock assessment developed by Froese et al. (2017) to the catch time series we extracted from the CONAPESCA database. Summarily, this method, which Froese et al. (2017) dub the “Simple CMSY,” applies a Bayesian model to fisheries-dependent time series to generate most-likely values for B/B_MSY over time. It does so based on the landings data combined with informed priors about the stock’s intrinsic population growth (r) and starting biomass (B₀ ). To better understand the theoretical and mathematical underpinnings of the Simple CMSY method, we recommend reviewing the original research by Froese et al. (2017) as well as a subsequent application of their methods to SSFs in Mexico executed by Giron-Nava et al. (2019). The resulting stock assessment profiles we generated using this method are included in Supplementary Materials - “ Supporting Dataframe 2 ”.

To test for a relationship between B/B_MSY and fuel intensity across the 39 “Stock-Years of Interest,” we employed three complementary statistical analyses. We reasoned that if there exists a strong (linear) relationship between B/B_MSY and fuel intensity across SSFs, it would be apparent via a simple correlation test and, likely, a simple linear model. Thus, we began with (i) a Pearson’s correlation test between estimated values for B/B_MSY and mean values for fuel intensity (Fuel Intensity_{mean s}). We then constructed (ii) a simple linear regression, where Fuel Intensity_mean values are predicted by B/B_MSY . Each Fuel Intensity_mean value is calculated as the arithmetic mean of all FIEs associated with that particular Stock-Year. Here, we compared B/B_MSY values to Fuel Intensity_mean values in an effort to control against the influence of uneven sample sizes across Stock-Years (sample sizes for each are listed in Table S2 ).

The third and final analysis we conducted was (iii) a mixed effects model where any individual estimate for fuel intensity is predicted by its corresponding value for B/B_MSY as a fixed effect, along with three random effects: gear type, genus, and year. One might reasonably expect that the fuel-use per kilogram of landed seafood on any given trip is borne from a number of factors beyond the fishery’s underlying stock biomass which could cause or covary with B/B_MSY . We know, for instance, that fishing patterns, fleet (over)capacity, management schemes, underlying ecosystem dynamics, weather patterns, and changes in climate are all important factors that can determine fishery outcomes (e.g., Worm et al., 2006; Pershing et al., 2015; Schuhbauer et al., 2017; Hilborn et al., 2020; Giron-Nava et al., 2021). For the most part, these dynamics lie beyond the scope of this paper, or the data are not available at this time. We can, however, begin control for some of this complexity — both observed and unobserved — by incorporating the additional information we have access to via the GCMP database about gear type, genus, and year (at the trip-level).

We controlled for gear type as a random effect based on findings from previous studies, which indicate that gear type plays an important role in the overall emissions borne from fishing activities (Parker and Tyedmers, 2015; Parker et al. 2018; Ferrer et al., 2021) . Likewise, we included genus and year as random effects since target taxa and unobserved changes over time likely play important roles in predicting fuel expenditures. We chose to model these three factors as random effects because, while we anticipate that they do matter, we are not explicitly interested in the fixed effects of gear type, genus, nor temporal factors on the observed discrepancies among individual estimates for fuel intensity. Rather, our primary interest lies in understanding the generalized relationship between B/B_MSY (a fixed effect) and fuel intensity (a continuous output) across gear types, target genera, and years.

We generated the mixed effects model using the “lme4” package in R (Bates et al., 2015), and utilized the R “lmerTest” package to generate associated p values (Kuznetsova et al., 2017). Notably, the “lmerTest” package applies “Satterthwaite’s method” (Fai and Cornelius, 1996) to estimate degrees of freedom, covariance, and significance values for unbalanced samples in mixed effects models. Thus, in contrast to our first and second analyses, it was not necessary to test B/B_MSY against Fuel Intensity_mean values here, due to the very nature of the mixed effects model itself. Instead, we compared estimates for B/B_MSY against individual estimates for fuel intensity, generated over 4,491 single-species fishing trips. (Data associated with each of these trips can be accessed in Supplementary Materials - “ Supporting Dataframe 3 .”)

We have included the results from our stock assessments in “ Supporting Dataframe 2 ”. Specifically, we report estimated values and confidence intervals for B/B_MSY and F/F_MSY for all 19 “Stocks of Interest” over time (2001–2019), as generated by the Simple CMSY method.

Stocks that are both overfished and experiencing overfishing are those with a B/B_MSY < 1 and a F/F_MSY > 1 (e.g., Giron-Nava et al., 2019; see also UW - Sustainable Fishing 101: https://sustainablefisheries-uw.org/seafood-101/overfished-overfishing-rebuilding-stocks/). Conversely, stocks that are fished “sustainably” are those that are neither overfished (B/B_MSY > 1) nor experiencing overfishing (F/F_MSY < 1). There are many cases in between, where stocks can be overfished but recovering (B/B_MSY < 1, F/F_MSY < 1), or not yet overfished but experiencing overfishing (B/B_MSY > 1, F/F_MSY > 1). With these definitions in mind, we found evidence of at least some overfishing in every year for which we were able to generate assessment data. We also found evidence that many of these stocks are currently overfished.

In 2019, the most recent year for which we were able to generate assessment data, we found that out of the 19 stocks, nine (47%) were overfished and experiencing overfishing (B/B_MSY < 1, F/F_MSY > 1); zero (0%) were overfished and recovering (B/B_MSY < 1, F/F_MSY < 1); two (11%) were not overfished but were experiencing overfishing (B/B_MSY > 1, F/F_MSY > 1); and eight stocks (42%) were being fished sustainably (B/B_MSY > 1, F/F_MSY < 1). Of the eight stocks that were fished sustainably that year, six of them were located in the LP, one in the LG, and one in the UG.

Our results for B/B_MSY , on which we’ve based our assessment of overfishing, are largely consistent with anecdotal evidence, ecological research (e.g., Lluch-Cota et al., 2007), and local ecological knowledge about the region’s fish populations (e.g., Sáenz-Arroyo et al., 2005). For example, as of 2016, Giron-Nava et al. (2019) found that in the GoC, 69% of the stocks analyzed (n = 121) were overfished and still being fished at unsustainable levels; 13% were overfished but recovering; 11% were not overfished but were being fished at unsustainable levels; and just 7% were fished “sustainably,” neither overfished nor experiencing overfishing.

A visual assessment of the data ( Figures 3A, B ) suggests that, the relationship between B/B_MSY and fuel intensity is inverse and exponential of some type, and we therefore report our results for B/B_MSY against log₁₀(Fuel Intensity) or log₁₀(Fuel Intensity_mean ) values where appropriate. We’ve elected to report the results of our first two analyses in terms of B/B_MSY versus Fuel Intensity_mean for reasons explained above (in Section 2.5), however, the interpretation of our results does not change if we compare values of B/B_MSY to demeaned values for fuel intensity ( Figure S1 ).

We found a statistically-significant negative correlation between B/B_MSY and log₁₀(Fuel Intensity_mean ) using a Pearson’s test for correlation, where r(df = 37) = -0.44 (p value < 0.01). Subsequently, we conducted a simple linear regression, where log₁₀(Fuel Intensity_mean ) was predicted by B/B_MSY with an R² _adj. = 0.17 (df = 37, p value < 0.01) ( Figure 3C ). Finally, the existence of an inverse relationship between B/B_MSY and fuel intensity was further corroborated by the results of our mixed effects model ( Figure 3D ), where B/B_MSY is a statistically-significant predictor of log₁₀(Fuel Intensity) (p value < 0.01), when accounting for gear type, genus, and year as random effects. According to this model parameterization, the association between B/B_MSY and log₁₀(Fuel Intensity) has a negative correlation of r = -0.36. Notably, B/B_MSY remained a significant predictor of both fuel intensity and Fuel Intensity_mean under a variety of controls and robustness checks ( Figure S1 ).