Relating a recruitment shift of Patagonian grenadier (Macruronus magellanicus Lönnberg) to large-scale environmental changes off Southern Chile

Abstract

The SE Pacific stock of Patagonian grenadier (Macruronus magellanicus) showed evidence of an abrupt reduction in recruitment after 2000. This drop exceeded expectations from changes in the spawning stock biomass (SSB), indicating a change in the stock-recruitment relationship (S-R). We evaluated whether variability in recruitment could be explained by concurrent changes in three environmental indices: sea-surface temperature anomaly (SSTA); southern oscillation index (SOI); and latitudinal position of the west wind drift bifurcation (WWDL). Continuous and discrete effects of these indices were tested as covariates in linear log-log and non-linear Ricker's S-R models. Discrete effects represented regime shifts detected in SSTA (1998), SOI (1998) and WWDL (1995). While SSTA was the only continuous variable found to be significant, the discrete 1998 regime shift supported the most informative model. The best Ricker's model considered a discrete intercept change in the same year: 1998. Although a spurious correlation between SSTA and S-R changes is possible, SSTA may be reflecting major physical or biological changes relevant to M. magellanicus juveniles in the SE Pacific.

Keywords:

• stock-recruitment
• west wind drift
• SSTA
• SOI
• climate change
• regime shift

Introduction

In fisheries science, fluctuations in marine fish abundance and resource productivity are explained by fishing mortality and recruitment variability, which is often related to survival in the early life history of fish. Several hypotheses have been postulated to explain recruitment variability around a mean that is expected to be driven by the spawning stock biomass (Myers & Barrowman 1996; Myers 1997). Most of these hypotheses are related to environmental variability and include the match-mismatch (Cushing 1975, 1990), starvation and predation (Hjort 1914; Lasker 1978, 1981; Rothschild & Osborn 1988; Purcell & Grover 1990), and advection (Parrish et al. 1983) hypotheses, as well as integrated frameworks such as the optimal environmental window (Cury & Roy 1989) and the triad hypothesis of Bakun (1996).

Assessing and incorporating environmental variability into stock-recruitment (S-R) models is expected to reduce uncertainty and improve S-R predictive power. Nonetheless, such variability has been incorporated in only a few S-R models actually used for management purposes (Walters & Collie 1988; Myers 1998; Francis et al. 2006). Most of these studies have focused on inter-annual fluctuations (‘recruitment variability’), while a more limited number have considered lower-frequency, high-amplitude recruitment shifts occurring at decadal or multi-decadal time-scales (Mueter et al. 2007; Wayte 2013).

Patagonian grenadier (Macruronus magellanicus Lönnberg) is a mesopelagic fish species that is widely distributed around the southern tip of South America (Fig. 1). Its main distribution range extends south of 29°S, ranging from 80°W in the Pacific to 56°W in the Atlantic Ocean (Arana 1970; Queirolo et al. 2008; Schuchert et al. 2010). Its main spawning grounds are located in the Pacific Ocean (44°10′S–46°20′S). Although the precise location and relative importance of nursery areas is poorly known, juveniles of ages 0+ and 1+ have been detected south of 40°S, in the coastal ocean, and within the extensive Patagonian Estuarine System (41°30′–47°S).

Figure 1 A, Geographic distribution range of Macruronus magellanicus around South America and main oceanic currents in the SE Pacific. B, Study area showing the 1×1° cells and specific geographical domains used to compute SSTA averages along the Chilean coast. WWD-B = west wind drift bifurcation, HC = Humboldt Current, CHC = Cape Horn Current, ACC = Antarctic Circumpolar Current, ONA = oceanic nursery area (hypothetical), ENA = estuarine nursery area, MSA = main spawning area.

Macruronus magellanicus is an important fishery resource that has endured substantial fishing pressure in Atlantic and Pacific waters. In the SE Pacific, very high exploitation rates were imposed by the pelagic purse-seine fleet that targeted juvenile aggregations off central Chile (36°S–42°S) between 1978 and 2003, averaging annual landings of 145,000 t (Cubillos et al. 2009). Currently, adult-oriented bottom trawl fisheries operate off central and southern Chile with average landings of 66,700 t. The 2010 Pacific stock assessment indicated a spawning stock biomass (SSB) below 20% of virgin conditions. Because fishing mortality has remained lower than reference limits during the last decade (Payá & Canales 2011), current overfishing and low SSB are probably related to a period of very low recruitment, starting in 2000 (Fig. 2). The recruitment shift occurring in 2000 could not necessarily be explained by the declining trend observed in SSB, rather we hypothesise that it represents a long-term recruitment shift, likely related to some persistent change in the environmental forces that affect M. magellanicus recruitment.

Figure 2 A, Time series of SSB. B, One-year-old fish recruitment of Macruronus magellanicus off southern Chile (1980–2010). Dotted lines are indicative of regime shifts identified in Table 1.

The objective of this paper was to evaluate the existence of a significant shift in M. magellanicus S-R during the last decades, and whether concurrent large-scale environmental changes could be inferred from observed anomalies in the sea-surface temperature (SSTA), the southern oscillation index (SOI) or the latitudinal position of the west wind drift (WWD) bifurcation. Links between SSTA, SOI and recruitment have been suggested in several marine species (Jacobson & MacCall 1995; Ottersen & Sundby 1995; Rutherford & Houde 1995; Fox et al. 2000; Galindo-Cortes et al. 2010), including species of the same genus such as hoki M. novaezelandiae in New Zealand (Livingston 2000; Bull & Livingston 2001). While the mechanisms explaining these links are unclear, it is probable that SSTA and SOI act as an index for a number of correlated factors affecting physical and biological processes at regional or ecosystem scales (Mueter et al. 2007).

The particular location of M. magellanicus main spawning grounds, just east of the WWD bifurcation (WWD-B), suggests that the latitudinal positioning of this bifurcation, where the Humboldt and Cape Horn currents are born, could affect the transport of eggs and larvae into suitable nursery habitats either in coastal or estuarine waters. Thus, although correlated to SSTA, the latitudinal position of the WWD-B (WWDL) becomes an alternative index that allows us to test some more specific hypotheses regarding the environmental forces that may affect recruitment in this population.

Materials and methods

Time series of spawning stock biomass and recruitment were obtained from the 2010 stock assessment model (Payá & Canales 2011). This is an integrated statistical catch-at-age analysis (Quinn & Deriso 1999) that assumes a single stock and population in the SE Pacific (Galleguillos et al. 1998; Machado-Schiaffino & Garcia-Vazquez 2011). We defined recruitment as the number of yearlings hind-cast at the beginning of the subsequent year (age 1 + ), and SSB as the total mass of spawning adults during the spawning season. Model parameters were estimated by maximum likelihood, contrasting model outputs with annual catches, fishery age composition, catch per unit effort, acoustic biomass index, and survey age-composition. Annual recruitment (R_t) was estimated through:

(1)

where is the average recruitment, and ε_t is the recruitment anomalies, with is the variance of recruitment and set equal to 0.6. A natural mortality of 0.35 yr⁻¹ was assumed constant across ages and years. Model data and AD Model builder template are available and can be solicited to the first author of this paper.

Sea-surface temperature anomalies (SSTA) were downloaded from the Integrated Global Ocean Services System (IGOSS) website (iridl.ldeo.columbia.edu). SSTA data (1981–2010) was averaged on a monthly basis across a 1×1° grid. We selected grid cells adjacent to the Chilean coast, between 36°S and 47°S (Fig. 1). This latitudinal range contains the main fishing, spawning and juvenile concentration grounds known for the species (Cubillos et al. 2009, Lillo et al. 2010). SSTA values corresponding to the post-spawning months (July–December) were averaged within years considering the overall distribution of the species north of 47°S and three specific spatial domains (Fig. 1), aimed to represent: 1. oceanic nursery areas (40°–42°S, hypothetical); 2. estuarine nursery areas (41°30′–47°S); and 3. main spawning area (44°–47°S). The latitudinal position of the WWD-B at 75°W (WWDL) was computed from isotherm distributions from average monthly SST following the procedures described by Gatica et al. (2009). The SOI was obtained from the US Climate Prediction Center (CPC, National Oceanic and Atmospheric Administration, www.cpc.noaa.gov/data/indices/sstoi.indices). Both WWDL and SOI data were subset and averaged to represent the July–December period within each year.

To reveal large-scale shifts in recruitment and environmental series, we used and compared results from two statistical approaches: 1. the chronological clustering method of Legendre and Legendre (1998); and 2. the Bayesian analysis for change points described by Barry and Hartigan (1993) and implemented in R by Erdman and Emerson (2007). Here, we conducted 5000 simulations, burning the first 50 samples from each chain and using priors of 0.2 for both the signal/noise ratio and the probability of detecting a shift.

To assess whether variability in S-R could be explained by selected environmental indices (EI), we considered two types of models. First, a family of linear models that related recruitment rate and EI, following the expression:

(2)

where a and b are regression parameters; and, second, a non-linear density-dependent model that followed the classical Ricker's (1954) model:

(3)

but allowed for changes in parameters a, b or both, conditional to the environmental shifts being tested.

Both continuous and discrete environmental effects were evaluated in the linear model, but only discrete effects were tested in Ricker's model to facilitate interpretation and avoid convergence problems. Discrete environmental effects were represented by categorical variables indicative of regime shifts being detected within each series. To model and account for temporal correlation, a mixed model approach was used for both linear and non-linear models (Pinheiro & Bates 2000). Errors were assumed to follow an auto-regressive structure of order one, which was the most informative among several tested ones. Model comparisons were based upon Akaike's (1973) information criterion (AIC).

Results

Both methods used to detect regime shifts yielded consistent results for all response variables with the exception of SSB. Here, the chronological clustering estimated a non-significant shift in 2002, while the Bayesian method estimated a 79% probability of a shift, except in 1982 (Table 1). This 1982 and 2002 interval bounds a period when SSB exhibited a consistent and declining trend (Fig. 2). SSTA exhibited a significant shift in 1998, regardless of the spatial domain being considered. This shift reflected the start of a cooling period that lasted for most of the remaining years in the series (Fig. 3). While no regime shifts were detected for SOI, a 2.1° poleward displacement of WWDL was detected in 1995, although it exhibited a posterior probability of only 32% (Table 1).

Figure 3 Annual averages of SSTA in coastal waters of the SE Pacific Ocean (35°S–47°S) from 1981 to 2010. Data intended to represent the first months of life of Macruronus magellanicus, from July to December. Colour scale on the right indicates temperature anomaly in °C.

Table 1 Detection of shifts in sea-surface temperature anomaly (SSTA), latitudinal position of the west wind drift bifurcation (WWDL), southern oscillation index (SOI), spawning stock biomass (SSB) and recruitment time series (1981–2010).

	Chronological clustering		Bayesian change point		Relative change
Response variable	Shift year	P value	Shift year	Posterior prob.	(Regime 2– regime 1)
SSTA
Fishing area	1998	<0.001	1998	0.89	−0.65 °C
Shelf nursery area	1998	<0.01	1998	0.68	−0.47 °C
Estuarine nursery area	1998	<0.01	1998	0.93	−0.79 °C
Spawning area	1998	<0.01	1998	0.92	−0.57 °C
WWDL	1995	<0.001	1995	0.32	+2.13 °S
SOI	1999	>0.1	1999	0.11	+0.88
loge(R)	2000	<0.05	2000	0.82	−0.87
loge(SSB)	2002	>0.1	1982	0.79	−0.78
loge(R/SSB)	1998	<0.05	1998	0.75	−0.25

Note: P values for chronological clustering correspond to the maximum likelihood ratio between autoregressive (order one) mixed models, fit with and without the corresponding shift. Shift years correspond to those that exhibited the maximum likelihood or the highest posterior probability to represent a shift within each series.

All linear models that related recruitment rate to SSTA averaged within different spatial domains showed a significant relationship (P < 0.05 for all models). Within this context, the most informative domains corresponded to the estuarine nursery (AIC = 25.4) and main fishing (AIC = 25.9) areas, and the least to the spawning (AIC = 26.6) and shelf nursery (AIC = 26.7) areas. No significant linear effects upon reproductive rates were found for either WWDL or SOI (P > 0.05 for both). Replacing continuous explanatory variables by a discrete 1998 regime shift (Fig. 4A) led to a significant improvement in the linear model (AIC = 23.7). The discrete effect of a 2006 regime shift (WWDL) was non-significant.

Figure 4 A, Observed relationships between hind-cast recruitment rates (R/SSB) of Macruronus magellanicus and SSTA in its main fishing area 1981–2008 and predicted values from an auto-regressive linear regression model. B, Relationship between hind-cast recruitment (R, in millions) and SSB (in tonnes) during the same period, and predictions from a non-linear conditional Ricker's (1954) S-R model. Year 2008 is highlighted in both panels as an anomalous year, although it was not excluded from any analysis.

Although the non-linear fit of Ricker's model was rather poor, the results were consistent with those obtained from the linear model. Here, the null hypothesis of constant parameters was the least informative one (AIC = 23.4), while the model that made parameter a conditional to the 1999 shift exhibited the highest performance (AIC = 16.7, Fig. 4B). The full model that estimated independent sets of parameters a and b before and after 1999 showed an intermediate performance (AIC = 20.6).

Discussion

This study presents evidence that the SE Pacific stock of M. magellanicus experienced a significant recruitment shift in 2000 which did not match any concurrent shift in SSB, but which exhibited a severe and progressive decline between 1996 and 2002. We found consistent evidence that the reproductive rate and the S-R experienced an abrupt change in 1998 (Fig. 4). This shift was concurrent with large-scale environmental changes reflected by SSTA in the same year (1998), particularly at the spawning grounds (44°S–47°S). Thus, we speculate that the 2000 drop in recruitment resulted from the combined effects of low SSB levels and adverse environmental conditions affecting early juveniles since 1998. Furthermore, anomalous R and R/SSB values observed in 2008 (Fig. 4) might represent a new shift that needs to be explored in the future.

Besides overfishing, the rare presence of juveniles north of 42°S during the last decades could be related to a limited access or reduced suitability of this area to juveniles. After 1995, the WWDL showed poleward displacement. This change could be related with favourable advection of early life stages into nursery areas, but it was not supported by our findings. Indeed, we found no significant relationships between WWDL and recruitment. We found, however, larger SSTA changes in the Patagonian Estuarine System and a tighter relationship between recruitment and SSTA variability here than in oceanic nursery areas. Therefore, it may be necessary to consider the relative importance of advective forces into both nursery areas, at scales and depths (0–100 m) that may be more relevant to early juveniles than the WWDL data available to us.

The concurrence and correlation we found in trends and shifts of SSTA and recruitment are not tests of causality, and may have resulted from simple synchrony between these variables. Testing these relationships with new datasets is a necessary step that has failed to yield significant results in many similar cases (Myers 1998). As a related example, none of the relationships initially found in M. novazelandiae between short-term variability in recruitment and climate (Livingston 2000; Bull & Livingston 2001) remained significant when another six years of data were added to the series (Francis et al. 2006). In the meantime, there is a need to identify which environmental factors, relevant to recruitment, are actually being indexed by SSTA and whether endogenous mechanisms may have induced disproportional and discrete recruitment shifts in S-R when a given threshold SSB was reached.

At least six regime shifts have been reported for the North Pacific since 1920, the last of them in 1999 (Breaker 2007). Although no evidence of a 1999 regime shift existed for the Humboldt Current System, available publications focused mainly on fluctuations in pelagic species (Chavez et al. 2003; Alheit & Niquen 2004). Our results, which include latitudes in the southern range of the Humboldt Current and beyond, suggest the need to revisit this issue to conduct a larger-scale analysis of SSTA and other time series in the SE Pacific area. Whether of natural or anthropogenic origin, this low frequency variability has the potential to induce large changes to both fisheries and marine ecosystems.

Acknowledgements

IFOP provided access to information, data and codes related to the stock assessment of M. magellanicus. LC and EN work was supported by the National Commission for Scientific and Technological Research (CONICYT), FONDECYT N° 1111006.