Environmental changes in the northern Gulf of California (NGC) are thought to negatively impact the primary productivity of the region, with adverse consequences on the development of eggs and larvae of numerous marine species (invertebrates and fishes, Álvarez-Borrego 1983, Silber 1990). Another factor considered to negatively affect marine populations is fishing activity by industrial shrimp and shark fleets, through retention of juveniles and pre-adults. This activity also impacts species of particular conservation interest in the region, such as totoaba (Totoaba macdonaldi) and marine vaquita (Phocoena sinus), due to their status as endangered species (Diario Oficial de la Federación de México 2010). The NGC and the Colorado River Delta currently constitute a natural reserve protected by the Mexican Federal Government and included since 1995 in the National Programme of Natural Protected Areas (Fig. 1).

The relationship between oceanic patterns and variability in fisheries has been documented previously (Hare and Mantua 2002). In the Gulf of California, fluctuations in sardine fishery are related to intradecadal events such as the El Niño Southern Oscillation (ENSO) (Lluch-Belda et al. 1986). Lluch-Belda et al. (1992) also observed that the Pacific Decadal Oscillation (PDO) index, has a significant effect on the sardine populations (Arreguín-Sánchez and Martínez-Aguilar 2004). Similar relationships occur for other species in the Gulf of California, including fish and shrimp (Lluch-Cota et al. 2007, Arreguín-Sánchez et al. 2017), primary producers (Espinosa-Carreón et al. 2004), totoaba (Lercari and Chávez 2007), giant squid (Nevárez Martínez et al. 2002) and pelicans (Sydeman et al. 2001). These reports suggest interactions between the physical environment and the ecosystem that affect species throughout their ontogeny (Lluch-Cota et al. 2007).

It also has been documented that large-scale oceanic climate indices such as the PDO index and the Multivariate ENSO Index (MEI) (Wolter 2018) have relevant effects on fisheries in the Gulf of California (Castro-Ortiz and Lluch-Belda 2008). The Pacific Ocean has experienced an interannual-to-interdecadal recurring pattern of ocean-atmosphere climate variability centred over the mid-latitude Pacific (Zhang et al. 1997, Newman et al. 2003, Yu et al. 2017), which has significant influences on the fishery (Litz et al. 2011, Newman et al. 2016). The PDO is characterized by the two phases of the sea surface temperature (SST) anomaly: positive (warm) and negative (cool) phases. During the positive phase, the water becomes cooler in the northwestern and central Pacific Ocean and warmer in the eastern Pacific Ocean. The SST anomaly pattern is reversed during the negative phase: warmer in the northwestern Pacific Ocean and cooler in the eastern Pacific Ocean (Mantua and Hare 2002). On other hand, the Pacific–North American (PNA) pattern affects coastal sea and continental surface air temperatures, as well as streamflow in major west coast river systems and terrestrial precipitation variability in North and Central America (Yu and Lin 2019).

The Gulf of California is a marine system with a complex environment. One of the greatest challenges for resource managers is to understand how much of the observed change in a population’s abundance can be explained by the environment. Understanding these relationships will enhance the probability of success in managing populations for both fishing and conservation. This is especially true if fishing and conservation efforts can be adapted based on observed variability and if species show specific responses to environmental factors in relation to life histories. Therefore, the focus of the present study is to estimate the role of environmental variables in explaining observed variability of the relative abundance of several commercial species in the NGC.

Catch per unit effort (CPUE) (kilogram/boat, k/b) was used as an indicator for the relative abundance of commercial species (Spare and Venema 1998). The time series of daily catch and fishing effort for commercial species were obtained from logbooks of the Regional Fisheries Offices (1995 to 2007). Species catches were listed by common commercial names which may include one or several species, such as chano (Chanos chanos, Micropogonias megalops), corvina (Cynoscion othonopterus, C. parvipinnis, C. xanthulus, C. reticulatus, C. nobilis), blue shrimp (Litopenaeus stylirostris), Pacific sierra (Scomberomorus sierra, S. concolor), ray (Rhinobatos productus, Dasyatis brevis, Gymnura marmorata, Urolophus maculatus, U. halleri, Myliobatis californica, Rinopthera steindachneri, Narcine entemedor), shark (Heterodontus francisci, Heterodontus mexicanus, Rhizoprinodon longurio, Mustelus spp., M. lunulatus, Sphyearna lewini, Squatina californica) and totoaba (Totoaba macdonaldi). For the totoaba, the catch time series included 1929 to 1971, when the fishery was active (Arvizu and Chávez 1972) and the fishing effort corresponded to a reconstructed time series (Lercari and Chávez 2007).

The selection of environmental variables affecting species considered in this study was based on information from scientific publications: time series of chlorophyll and primary productivity (Espinosa-Carreón et al. 2004), SST (Aragón-Noriega 2007), ENSO (Lluch-Cota et al. 2007, 2010), the Colorado River discharge (Rowell et al. 2005, Pérez-Arvizu et al. 2009, Aragón-Noriega y Calderón-Aguilera 2000), precipitation (Valdez-Muñoz et al. 2010), the PDO index and the MEI (Lluch-Cota et al. 2007, 2010) and the PNA pattern, (Cabello-Pasini et al. 2003), and the sources are shown in Appendix 1. Generalized additive models (GAMs) were used to identify variables that explain changes in CPUE because they allowed several environmental variables acting on stock abundances (CPUE) to be incorporated during a single time. The resulting GAMs for each group were chosen considering the lower AIC value (Akaike 1973), and were also validated by the generalized cross validation estimator (GCV), used as a measure of goodness of fit period (Wood and Augustin 2002, Venables and Dichmont 2004). The GAM model was based on R script (R Core Team 2018) version 3.5.2, using the mgcv package version 1.7-5.

Broad-scale climate indices not only affect weather and climate variability worldwide but have also been demonstrated to considerably influence ecological processes (Stenseth et al. 2003). To identify possible biological effects related more strongly to long-scale climate patterns than seasonal factors, and to avoid masking biotic response by seasonal effects, CPUE time series were transformed using natural logarithm (x+1), and were also standardized by z-transformation to facilitate visualization (Blanchard et al. 2010). Then we conducted a seasonal and trend decomposition with repeated locally weighted scatterplot smoothing (LOESS), using the stl function in the R package (Zuur et al. 2007). We used a large span width of periodic signal (blue shrimp, chano and sierra, 6 months; sharks and corvina, 5 months; rays, 8 months; totoaba, 2 years) to provide an overview of long-term trends. We also used the trend component and verified for uncorrelated noise processes.

De-seasonalized significant environmental and climate time series derived from GAM models were used to test cross-correlation analysis between them and CPUE of all commercial groups, using the ccf function in R (Walton et al. 2017). This function is an estimate of the association between two variables (CPUE and the physical environmental variables) at the significant temporal leads and lags (George et al. 2005). Temporal lags in the relationship between predictors and the response variable were accounted for, because it has been reported that climate indices influence recruitment success of fish species several months and even years later (Kashkooli et al. 2017). Positive time lags represent a positive correlation between driver and response and could indicate causality. Negative correlation would be represented by time lags opposed in the driver and response variables (Walton et al. 2017).

For the predictive model, seven groups of exploited species were considered, including nine fish species (including the totoaba), sixteen elasmobranchs and one crustacean species. Table 1 shows model parameters and the significant environmental variables that corresponded to deviances in abundance for different species. Figure 2 shows the observed CPUE values of the species analysed, as well as the resulting fittings through the GAM for each population. The GAM coefficients and the respective standard errors for all populations analysed in this study are presented in Appendix 2. According to the results obtained, SST was statistically significant with respect to the variability in CPUE, except for the totoaba, for which the Colorado River discharge (CRD) explained the greatest variability in abundance. Other environmental variables that showed significant effects on several species were the PDO and the CRD.

For both sharks and rays the statistically significant environmental variable that explained most of the deviance in CPUE was SST, with 23.69% explained for sharks and 52.10% for rays (Table 1). Based on these results, it was further observed that the highest values of relative abundance (based on CPUE) occurred within an SST range of 22°C to 29.5°C for sharks and 16°C to 20°C for rays (Fig. 3A and B, respectively). The PDO index explained 18.36% of the deviance for sharks and 5.10% for rays. For both, greater abundance was observed during the cold phases of the PDO (Fig. 3C and D, respectively). The higher abundances for both sharks and rays occurred between the cold to neutral phases (Fig. 4A and B, respectively). The total deviance explained using the GAM model is given in Table 1, while Figure 2A and B show observed vs estimated CPUE values.

For blue shrimp, the environmental variable with the most significant impact was SST, which contributed 32.7% of the explained deviance, followed by the MEI, which contributed 12.4%. For SST, the greatest abundance of shrimp was observed within the range of 24°C and 30°C; for the MEI, higher CPUE values were observed when the anomaly remains within a range of ±1°C (Fig. 3E and F, respectively). During El Niño event periods, low CPUE values were recorded. The total deviance explained by the GAM model was 45.10% (Table 1). Figure 2C shows the observed CPUE values vs. the CPUE values estimated by the GAM method, and Figure 4C shows the standardized values of the CPUE and the significant environmental variables related to the blue shrimp.

The most significant environmental variable affecting the Pacific sierra was SST, contributing 67.10% of the explained deviance for variability in CPUE, followed by the PNA index, contributing 7.7% (Table 1). For SST, the highest abundance values were observed within the range of 22°C to 27°C (Fig. 3G), coinciding with the negative phases of PNA pattern (Fig. 3H). The GAM model explained a total deviance of 74.80% of the variability in serranid abundance, including both SST and the PNA pattern (Fig. 2D). Figure 4D shows the standardized values of CPUE and the significant environmental variables related to the Pacific sierra.

For sciaenids, including both chanos and corvina, SST explained most of the variability in CPUE. SST explained 63.10% of the deviance in CPUE for chanos and 67.30% for corvina. The results further indicate a positive correlation between SST (up to 25°C) and changes in the abundance of chanos and a negative correlation for SST (values >21°C). Maximum abundance occurred in the range of SST from 19°C to 24°C (Fig. 3I). For the corvina, a negative correlation between abundance and SST was identified for some ranges. CPUE of the corvina showed no changes at temperatures below 18.5°C. When SST was greater than 19°C, the relationship between this and the CPUE of the corvina showed an inverse relationship (Fig. 3J). The CRD was another environmental variable that significantly affected the chanos and contributed 13.2% to the explained deviance of CPUE variability (Fig. 3K); for the corvina, the CRD contributed 13.9%. The CRD showed a positive correlation for the chano, but for the corvina, for discharges lower than 9 m³ s^–1 of water of Colorado River, CPUE of the corvina showed an inverse relationship with the CRD; for discharges higher than 9 m³ s^–1, the abundance of the corvina showed a direct relationship with the CRD (Fig. 3L). The PDO index was also statistically significant for chanos, with 7.7% of the explained deviance (Fig. 3M), while PDO contributed 5.9% of the explained deviance for corvina (Fig. 3N). In both cases, the negative phase of the PDO resulted in favourable conditions for the abundance of both chanos and corvinas. Figure 4E and F illustrate the standardized values of CPUE and the environmental variables for the chanos and corvina, respectively. Finally, the GAMs resulted in a total explained deviance of 84% for chanos and 87.10% for corvina (Table 1, Fig. 2E and F).

For totoaba, the GAM suggests that the CRD explained most of the variability in CPUE, with 42.0% of the deviance, followed by the PDO with 15.2%. The total explained variation in CPUE was 57.20% (Table 1, Fig. 2G). The CRD exhibited a significant and positive relationship with abundance of totoaba at flows lower than 400 m³ s^–1. The highest abundances were observed when the CRD flows ranged between 200 and 500 m³ s^–1 (Fig. 3O). The PDO showed a trend that was similar to that of CPUE (Fig. 3P). These results suggest that the totoaba stock responds positively to the warm phases of the PDO. Figure 4G shows the standardized observed values of CPUE and the significant environmental variables related to GAM for totoaba.

We cross-correlated the CPUE of elasmobranchs, shrimp, sierra and sciaenids from the NGC with SST, the CRD and climate indices. Cross-correlations (Table 2) indicated that elasmobranch CPUE significantly decreased with SST with no time lag. Shark CPUE was positively correlated with SST anomalies three years earlier (lag of –3 year). Shrimp CPUE was also negatively correlated with SST anomalies and with no time lag, while the Pacific sierra was positively correlated with a time lag of –2 years. (Table 2). Abundance of chano and corvina were negatively correlated with SST, with a three- and one-year time lag, respectively (Table 2).

The long-scale climate PDO index was negatively correlated with elasmobranch CPUE, but with no time lag for sharks and a one-year time lag for rays (Table 2). The PDO was significant for sciaenids: negatively for chano and corvina, and positively for totoaba. Cross-correlation shows a two-year time lag for corvina and no time lag for chano. For totoaba the positive correlation also indicates a two-year time lag with the PDO (Table 2). Blue shrimp was the only commercial species that showed a significant correlation with the MEI. It was positive but showed a five-year time lag, while for the PNA pattern the Pacific sierra was the only one that significantly correlated, being positive and with two-year time lag (Table 2). Finally, the Colorado River discharge was positively correlated with all sciaenids, but with different time lags: 1 year for chano, –1 year for corvina and no time lag for totoaba. (Table 2).

According to our results, six of the seven groups under study demonstrated a high correlation between abundance and SST, as shown by several authors for sharks in the Mexican Pacific (Soriano Velázquez et al. 2006). This correlation support reports indicating that the SST plays an important role in seasonal migration for several elasmobranch species in coastal waters (Talent 1985, Wallman and Bennett 2006).

For rays, fisheries in the NGC have increased in recent years, but little research has been published regarding seasonal changes in abundance related to temperature (Cudney-Bueno and Turk-Boyer 1998). In both cases, considering the species mobility capacities, the literature suggests that the relationships with SST are related to habitat preferences.

For the blue shrimp, SST is the most significant environmental variable affecting stock abundance in the NGC, coinciding with fishing areas and greater abundances, and particularly with the reproductive process (Aragón-Noriega 2007).

For the above three study cases, cross-correlations suggest that temperature effects occurred at seasonal level (higher cross-correlations at time t=0), probably associated with ENSO effects and with mesoscale variability (McClatchie 2014).

With respect to the Pacific sierra (Scomberomorus sierra), Medina-Gómez (2006) and Valdovinos-Jacobo (2006) suggested that migratory behaviour related to feeding and reproduction is governed by changes in temperature, which favoured fishing activities (from October to May). Higher Pacific sierra CPUE coincided with the minimum average values of SST (colder conditions), and the cross-correlation of –2 years with SST is probably related to recruitment success, as shown in the fishery. In contrast, the behaviour of S. concolor, which is an endemic species in the NGC, is governed by the dynamics of the Colorado River Delta, particularly in relation to the reproductive process (Valdovinos-Jacobo 2006), suggesting a synchrony with warm waters favouring larval development (Moser et al. 1974).

It is well known that sciaenid fish such as chano (milkfish) and corvina (croaker) inhabit coastal waters, estuaries, rivers and deltas mainly during the reproductive season (Chao 1995). In the NGC, the milkfish shows spawning reproductive aggregations coinciding with warmer waters (CONANP 2007), as shown in Figure 3I, while the lag observed by cross-correlation suggests that SST is related to adults during the reproductive process (Table 2). In contrast, corvina (Cynoscion othonopterus), which is widely distributed in the Gulf of California, shows an inverse relationship with SST (Fig. 3J), coinciding with that reported by Chao (1995), while cross-correlation indicates higher value with only a one-year time lag (Table 2), probably reflecting the wide distributions of the species and habitat preference during the reproductive migration period in spring (Pérez-Valencia et al. 2012).

The GAM model indicated that shark and ray catches in the NGC significantly decreased in the same year as a high PDO (and increased after a low PDO anomaly) (Fig. 3C and D respectively). The influence of the PDO index during the shark fishery season could be explained by the fact that SST decreases because of the influence of cold weather and strong winter winds in the northern Pacific Ocean (Castro-Ortiz and Lluch-Belda 2008). Differences between the high CPUE for ray in spring and summer found by Ramirez-Amaro et al. (2013) and our results suggest that considerable fishing effort may be opportunistically directed.

The influence of the PDO in the dynamics of the NGC has been described by several authors (e.g. Hare and Mantua 2002), who pointed out that the PDO’s positive phase is characterized by a higher SST in the Alaskan current in the summer months. These events increased the CRD (through the drainage of dams) during the periods 1979-1981 and 1983-1987 (Lavín and Sánchez 1999, Rowell et al. 2005) (Fig. 5), and were associated with the increase in abundance of marine species in the NGC, including the corvina, whose abundance has been increasing since 1992 (Rowell et al. 2005). GAM models for chano and corvina show a negative correlation between CPUE and PDO anomalies (Fig. 3M and N).

Totoaba stock responds directly to the warm and cold phases of the PDO, showing a significant delay with PDO anomalies of –2 years. (Table 2). In this regards, Fischer et al. (1995) indicate that wind patterns and upwellings affect reproduction (Fig. 3P). According to Cury and Roy (1989), the synchrony between reproductive success and upwelling intensity can promote the success or failure of fish larva survival. In general, totoaba CPUE coincided with large-scale annual PDO and CRD trends (Fig. 4G).

Castro-Ortiz and Lluch-Belda (2008) reported the influence of the PDO on the declining of blue shrimp catch between 1988 and 1991, related to the cold period (Lavín and Sánchez 1999), and a steadily rising trend that peaked in 1996-1997. Our results suggest that longer blue shrimp reproductive periods occur when MEI anomalies are within the range of ±1 (called the neutral phase, Figs 3F and 4C), decreasing with MEI anomalies greater than ±1. Though blue shrimp CPUE is positively correlated with the MEI with a time lag of –5 year (Table 2), there is no evidence of the process behind this correlation.

High catches of sierra coincide with the negative anomaly of the PNA pattern, associated with colder conditions. It has been suggested that the effect of the PNA pattern occurs through changes in the displacement of air masses. and therefore changes in temperature and precipitation (Wallace and Gutzler 1981), and is frequently associated with El Niño-La Niña events (Wise 2012). The GAM model suggests that higher stock abundance values occurred during the negative and neutral phases of the PNA pattern, as evidenced by adult fishes having a greater presence in cold temperatures (Fig. 3H).

Cudney-Bueno and Turk-Boyer (1998) reported that corvina abundance seemed to increase with CRD increments during the period 1979-1988. This increase also corresponded to a decreasing abundance of totoaba, suggesting a species replacement process because of their similar feeding and reproductive niches. Our results in GAM models prove a positive correlation between sciaenid CPUE and CRD (Fig. 3K, L and O respectively), but with different time lags between species (Table 2). Correlations found between CRD and CPUE in sciaenids in the NGC are explained by the fact that chano, corvina and totoaba are endemic species that use the delta region for reproductive and breeding aggregations because they require estuarine conditions. The importance of CRD as an essentially environmental factor have been well documented (Gillanders and Kingsford 2002), while Rowell et al. (2005, 2008) suggested that higher CRD would increase spawning and nursery habitats, thus benefitting recruitment. Differences in years’ time lags could be explained as particular adaptive processes according to life histories.

Our analysis of totoaba only cover the period in which totoaba fishing was permitted (1925 to 1975). During this period, the CRD represented an environmental variable with a strong and significant effect on stock abundance. At present, the totoaba is considered a critically endangered species by the IUCN red list (Findley 2010), and the probable causes mentioned are changes in the CRD, spawning and nursery habitat losses, commercial overfishing and illegal captures, all of them interfering with the success of recovery management measures (Cisneros-Mata et al. 1995, Rowell et al. 2008). Since 1943, the CRD has been considerably reduced, coinciding with a very marked decrease in the abundance of totoaba. Unfortunately, present stock regulations (permanent closure) do not permit further investigation into the current effects of climate variables (due to the absence of abundance data), since it is expected that the habitat has changed over the last 30 to 40 years.

The NGC is a highly dynamic ecosystem influenced by several environmental variables that affect species differently in accordance with their life history. However, with the exception of the totoaba, temperature appears as a key environmental variable affecting stock abundance. PDO influences elasmobranches and sciaenid fishes, suggesting effects on species with high mobility, while the CRD affects sciaenid fishes, suggesting, in accordance with the literature, an impact on the reproductive habitat and/or behaviour. Because of its status as a critically endangered species, results are especially relevant for the totoaba. In general terms, the time lags observed between environmental variables and stock abundance seem to be explained through the traits or processes of the life histories of the species and their adaptation to habitat conditions. Our results demonstrate that environmental variables can explain some changes in stock abundance and must be taken into account in the expectations of the management policies of exploited stocks, particularly those aimed at stock recovery of endangered species.