Skillful prediction of tropical Pacific fisheries provided by Atlantic Niños

Tropical Pacific upwelling-dependent ecosystems are the most productive and variable worldwide, mainly due to the influence of El Niño Southern Oscillation (ENSO). ENSO can be forecasted seasons ahead thanks to assorted climate precursors (local-Pacific processes, pantropical interactions). However, owing to observational data scarcity and bias-related issues in earth system models, little is known about the importance of these precursors for marine ecosystem prediction. With recently released reanalysis-nudged global marine ecosystem simulations, these constraints can be sidestepped, allowing full examination of tropical Pacific ecosystem predictability. By complementing historical fishing records with marine ecosystem model data, we show herein that equatorial Atlantic Sea Surface Temperatures (SSTs) constitute a superlative predictability source for tropical Pacific marine yields, which can be forecasted over large-scale areas up to 2 years in advance. A detailed physical-biological mechanism is proposed whereby Atlantic SSTs modulate upwelling of nutrient-rich waters in the tropical Pacific, leading to a bottom-up propagation of the climate-related signal across the marine food web. Our results represent historical and near-future climate conditions and provide a useful springboard for implementing a marine ecosystem prediction system in the tropical Pacific.


Introduction
Before COVID-19, global fisheries supplied society with ∼96 million tons of fish. With 39 million people engaged in the primary sector and 4.6 million fishing vessels in operation, the fisheries industry exported USD 164 billion annually [1].
Despite the enormous anthropogenic pressure these numbers reveal, year-to-year variations in marine resources from certain oceanic regions are still largely determined by Earth's natural variability [2,3]. This is the case of the tropical Pacific, one of the most biodiverse and productive marine areas worldwide, where El Niño Southern Oscillation (ENSO) exerts a paramount influence [4][5][6]. During El Niño events, weakened trade winds allow the penetration of warm waters over the central and eastern tropical Pacific, leading to suppressed equatorial and coastal upwelling of nutrient-rich waters. During La Niña these conditions reverse, inducing colder sea surface temperatures (SSTs), enhanced near-surface nutrient supply and a subsequent response of marine life [7][8][9][10][11].
Boreal spring NTA and autumn Indian ocean SST forcings affect the tropical Pacific mainly through oceanic temperature advection [19,29,30]. Although for both drivers eastward-propagating oceanic Kelvin waves are excited, in the first case SST and thermocline depth anomalies are triggered over the central Pacific [19], and in the second the main effect of Kelvin waves is through constructive interference leading to temperature advection processes [30]. Distinctively, equatorial Atlantic boreal summer anomalous SSTs (a.k.a., Atlantic Niños/Niñas) produce changes in the Pacific equatorial winds that excite eastward-propagating upwelling oceanic Kelvin waves. These waves effectively reinforce anomalies in the eastern Pacific, thereby altering the thermocline and the upwelling of ocean deep waters [31,32]. The settled environmental conditions are, in all cases, conducive to ENSO events in the ensuing seasons, but their impact on equatorial and coastal upwelling over the eastern Tropical Pacific, where highly-productive upwelling-dependent marine ecosystems are located [1], is different [15,23]. Thus, it is important to assess their relevance for marine ecosystem prediction, a fracture that has remained uncertain for several reasons.
Scarce, short or too local historical biological records often hinder the identification of robust statistical links between climate predictors and yield responses [33], especially when considering vast areas. Dynamical prediction by means of effectively coupled climate, biogeochemical and ecosystem models has been proven successful over multiple coastal and open sea regions [8,[34][35][36][37]. In the tropical Pacific, Earth System Model (ESM) simulations forced by observed climate have shown that multiyear predictability of tropical Pacific marine primary productivity is attainable [8,36,38]. These forced simulations have the benefit of reducing large SST biases present in free (unforced) runs [15,[39][40][41], thus allowing identification of local and remote predictability sources for eastern tropical Pacific marine ecosystems (e.g. Atlantic Niños).
With recently released global marine ecosystem models forced by observed historical climate from the Fisheries and Marine Ecosystem Model Intercomparison Project (FishMIP [42]), rigorous model-based statistical prediction can be implemented by considering the entire marine food web and its associated fisheries. Released fished and non-fished Fish-MIP historical simulations not only permit the effective isolation of climate-driven ecosystem responses, but also characterize the top-down control of this anthropogenic forcing.
Given the above, we evaluate herein the potential for prediction of tropical Pacific marine ecosystems and fisheries by assessing the precursory role of global SSTs. To this end, we utilize atmospheric reanalysis data, historical catch records and FishMIP simulations encompassing global marine ecosystems.

Observational and historical ESM simulation datasets
Annual catch data from the Food and Agriculture Organization of the United Nations (FAO) major fishing areas Eastern Central Pacific (#77) and Southeastern Pacific (#87) for the period 1950-2014 were retrieved from the Sea Around Us Project (SAUP) database (www.seaaroundus.org/data/#/fao). These are a combination of official reported data from FAO FishStat database (www.fao.org/fishery/statistics/en) and reconstructed estimates of unreported data [43]. All types of marine species and fishing methods were considered altogether in the SAUP timeseries, including catches from exclusive economic zones. To retain interannual variability and minimize the potential impact of external factors (e.g. fishing effort, climate change impact) on historical catch data, a high-pass Butterworth filter with 11 year cut-off period was applied to the time series.
Observational 1 • longitude-latitude SST data were retrieved from the Hadley Centre Sea Ice and Sea Surface Temperature data set (HadISST [44]). Simulated physical and biogeochemical data for the historical period  were obtained from the Geophysical Fluid Dynamics Laboratory (GFDL) Carbon Ocean Biogeochemistry and Lower Trophics run (COBALT [45]). Physical-biological parameters from this simulation (so-called GFDL ocean reanalysis within the FishMIP community) and additional datasets are fully described in table 1. The COBALT ocean-ice simulation, performed with the Modular Ocean Model v4p1 [46], has a 1 • horizontal resolution (1/3 • along the equator) and 50 vertical layers with 10 m increase in the top 200 m. The simulation is forced by air-sea fluxes from the CORE-II data set [47], a blend of the National Centers for Environmental Prediction (NCEP) atmospheric reanalysis [48] and satellite data. Whereas sea surface salinity is restored to observations in GFDL-COBALT due to drift issues, SST and biogeochemical outputs are free running (hindcasts [45,49,50]). Full information of the GFDL-COBALT validation against physicalbiological observations (e.g. mixed layer depth, nutrient concentration, phytoplankton/zooplankton biomass, carbon/energy fluxes) is available in Stock et al [45].
In line with CORE-II forcing on GFDL-COBALT, NCEP reanalysis was selected to analyze atmospheric
EcoOcean is a complex spatial-temporal global ecosystem model that comprehensively simulates interactions within all trophic levels and species groups using the Ecospace [52] spatial-temporal model at its core. The depth dimension is implicitly considered in the model by preferent habitat patterns and food-web interactions (one dimensional). EcoOcean v1 uses SAUP effort data [53] spatially distributed across large marine ecosystem regions as fishing forcing. Effort data combine fleet information (e.g. number and type of boats) and dedicated time for fishing. The forcing is applied as a mortality rate for functional groups in the simulation through the removal of existing biomass per unit of time. The model is able to replicate historical catch successfully on a global scale. A full description of EcoOcean is available in Christensen et al [54].
BOATS is a global ecological model that determines size spectra of fish biomass as a function of net primary production and local temperature. Biomass production in higher trophic levels is limited by available photosynthetic energy, temperature-dependent growth rates and trophic efficiency scaling. The model enables a realistic representation of biological and ecological processes, allowing evaluation of the links between climate, biogeochemistry and upper trophic dynamics. BOATS employs a bioeconomic approach based on SAUP catch price data to determine spatial-temporal changes in fishing efforts. The vertical dimension is a single surface-integrated layer and movement of fish is disregarded. Full information on the model architecture is available in Carozza et al [55].
Macroecological is a static-equilibrium model that relies on predator-prey mass ratios, transfer efficiency and metabolic demands dependent on body mass and temperature to predict size and distribution of marine consumers (as a whole) at any given time and location. It has a single vertical (surfaceintegrated) layer and movement of fish is disregarded. Extended information on the model is available in Jennings and Collingridge [56].
Seven different simulations were analyzed in this study. These include fished and unfished runs with and without diazotrophs dynamics (table 2). All historical runs  are forced by GFDL-COBALT simulation, whereas an EcoOcean 2021-2054 projection is driven by IPSL-CM5A-LR RCP8.5 [57]. In the latter, fishing forcing is kept constant at 2005 rates [42]. Although GFDL-ESM2M was also available as driving forcing for 2021-2054, it was eventually disregarded due to its inability to simulate a connection reminiscent of the observed Atlantic-Pacific teleconnection subject of this study [58]. Complete information on simulation settings, acronyms and data sources are available in table 2. All common FishMIP model output variables considered are detailed in table 1 and are publicly available through the FishMIP global model repository [51]. The most detailed model simulation (EcoOcean) forced by historical fisheries effort and which includes diazotrophs dynamics is used as reference throughout the manuscript (hereafter EcoG-Fis), as it provides the closest possible approximation to reality.

Physical-biological indices and statistical methods
Monthly anomalies of physical-biological variables were calculated by removing their total period average to the time series of each grid point. Seasonal (4 month) anomalies were computed as aggregated monthly anomalies, with linear trends effectively removed to minimize the impact of anthropogenic climate change [59].
Physical-biological indices were determined by spatially averaging seasonal anomalies over regions of interest, e.g. Niño3 (150 . All indices were normalized by dividing the time series by their standard deviation. Whereas Atl3 was computed for SST on boreal summer (July-September; JJAS) to represent Atlantic Niño variability during its peak intensity [22] (a.k.a. Atlantic Equatorial Mode), Niño3 and global seasonal anomalies of physical-biological variables (cf table 1) were calculated as rolling 4 month averages between monthly lags −8 to 24, where lag 0 corresponds to JJAS, lag 1 (−1) to JASO (MJJA) of the same calendar year and so forth (back).
Statistical significance for correlation outputs (figures 1, 3 and 4; 95% confidence interval) was evaluated using a Student-t test that accounts for the autocorrelation of the time series and the calculation of effective degrees of freedom [60].

S4CAST model
The SST-based Statistical Seasonal forecast model (S4CAST) is a novel statistical seasonal model that utilizes SST data (predictor) to forecast the response of any given impact variable (predictand). To this end, the model performs a maximum covariance analysis (MCA [61]) between a predictor Y field (SST anomalies) and a predictand Z field for any given region, time period and lag.
MCA is based on single value decomposition applied to the non-square covariance matrix (C) of both Y and Z fields. The method calculates linear combinations of the time series of Y and Z (a.k.a. expansion coefficients U and V, respectively) to maximize C.
A leave-one-out cross-validated hindcast based on MCA outputs is then performed [62]. For this purpose, predictor/predictand data combinations from all years except the year being predicted at each training stage are selected for model training and validation [62,63]. A Monte Carlo test with 1000 random permutations is applied to MCA (U vs. anomalous field correlations) and cross-validation of skill scores for hypothesis testing (95% confidence interval). The Pearson correlation coefficient is selected as a qualitative indicator for model skill (sign of response).
Further details on the S4CAST statistical model framework are available in Suárez-Moreno and Rodríguez-Fonseca [64].

Global SST relationship with tropical Pacific fisheries
We find year-to-year variations of historical global SSTs and eastern Pacific total catch records from the SAUP database (FAO areas 77 and 87) robustly correlated for the period 1950-2014. The years when observed fish catches are higher than usual synchronize with well-developed La Niña-like SST anomalies over the tropical Pacific and elsewhere [11,65] ( figure 1(A)). Lagged SST anomalies from the previous calendar year depict a preceding Atlantic Niño signal, together with a developing La Niña over the Pacific ( figure 1(B)).
Likewise, we consider EcoOcean model historical simulation in a fishing forcing scenario (EcoG-Fis) for the same analysis. We first characterize the tropical Pacific in terms of biomass density of commercial species. As expected, this is the most productive and fluctuating region worldwide in the model (figures 1(C) and (D)). Despite the slightly smaller area and shorter period considered (FishMIP coverage; 1971-2004), results very similar to observations are obtained from the model (figures 1(C) and (D)). Apart from equatorial Atlantic SST anomalies and intrinsic 2-7 years Niña-Niño intermittency, no other significant relations with tropical SST anomalies (e.g. NTA, Indian) are found, either in observations or in simulations during the years preceding annual-averaged total catch anomalies (figure 1 and supplementary material figure S1 (available online at stacks.iop.org/ERL/16/054066/mmedia)).
Correlation between the Atlantic Niño index (Atl3 JJAS SST; hereafter Atl3-Had) and the tropical Pacific total catch corroborates the lagged one year Atlantic-Pacific relationship in the model and observations (supplementary material figure S2). As expected, this relationship substantially weakens and changes sign in subsequent year lags. In agreement with Rodríguez-Fonseca et al [16], a slightly stronger Atlantic-Pacific connection is found for 1971-2004 compared to 1950-2014.

Statistical prediction of tropical Pacific total catch using equatorial Atlantic SSTs
Owing to constraints in historical fishing records (averaged annual data over vast areas), we focus on gridded 1-degree monthly FishMIP modeling outputs for fisheries prediction. Again, we consider the most detailed model forced by historical fisheries efforts (EcoG-Fis), which provides the closest approximation to reality (table 2). We then use the Sea Surface temperature-based Statistical Seasonal forecast model (S4CAST [64]) to perform a cross-validated hindcast of simulated tropical Pacific marine yields: total catch at sea (tc) and total system carbon biomass density (tsb).
For all monthly lags, Atlantic Niño emerges as the main predictor of tropical Pacific total catch anomalies, which can be hindcasted over large-scale areas and up to three years in advance (figures 2(A)-(E) and supplementary material figure S3). The quantification of the surface area of the domain with skillful prediction returns values of up to 3.4 million km 2 at lag 20 (figure 2(F)). If total system carbon biomass density is considered as predictand (Z) instead of total catch, final results are very similar but shifted ∼9-12 months back in time (compare supplementary material figures S3 and S4), with a maximum of 4.2 million km 2 of surface area with skillful prediction at lag 8 (figure 2(F)).

Atlantic-Pacific physical-biological mechanism
To establish a physical-biological mechanism explaining the results above, we evaluate tropical Atlantic (Atl3) SST impact on Pacific (Niño3) region for different area-averaged physical-biological variables and monthly lags. The simulation selected for this purpose is again EcoG-Fis, although other marine ecosystem models (BOATS, Macroecological) and scenarios are also considered (cf tables 1, 2 and Tittensor et al [42]).
First, we check whether the Atlantic-Pacific physical teleconnection is realistically represented in GFDL-COBALT. This is motivated by the fact that SSTs are not restored to observations in the simulation and equatorial wind fields in CORE-II depict inaccuracies [50,67]. Particularly, inconsistencies in wind speed, humidity, and solar radiation have been reported between NCEP and CORE-II [47]. In spite of this, the good performance of GFDL-COBALT representing Atlantic Niño variability and associated impacts is confirmed in figure 3(A) and supplementary material figure S5, where results obtained are very similar to observational studies [16,28]. For a warming in the equatorial Atlantic, a strengthening of the Pacific Walker Circulation is observed, together with increased subsidence and anomalous easterlies over the western equatorial Pacific, features that promote subsequent La Niña conditions (supplementary material figure S5). Similarly, previous studies have shown that anomalous easterly winds over the western equatorial Pacific promote the propagation of equatorial upwelling Kelvin waves [23,31,32,68], thus favoring the emergence of nutrients from the seabed.
The linear correlation value between Atlantic Niño indices (Atl3 JJAS 1971-2004) in GFDL-COBALT and HadISST is 0.87, whereas that for Niño3 (NDJF) is 0.96 (p-values < 0.05). The fact that both SST variabilities are significantly governed by the movement of water in response to wind variations may be behind the good agreement between GFDL-COBALT and HadISST [39,49].
Lagged-correlations for biological variables in figure 3(A) enable characterization of the bottom-up propagation of Atl3-related signal (hereafter Atl3-GF) across the Niño3 region food web. Closely synchronous to the growth of cold nutrient-rich SST anomalies ( figure 3(A)), large phytoplankton production (lphy), total system (tsb) and total consumer (tcb) carbon biomass densities grow rapidly before monthly lag 0, and decay thereafter at different rates. The different declining rates are due to distinct variable definitions and growth rates of species included in each category. While lphy is large-phytoplankton dependent, tcb is dominated by zooplankton and tsb is a blend of both (cf table 1). Consistent with previous studies [69], we find no SST response by small phytoplankton (not shown). Variables negative lag behavior can be partially explained by construction of Atl3 index at its peak season instead of its growing stage (supplementary material figure S6), and the fact that ocean chlorophyll tends to slightly precede ENSO-related SST anomalies due to wind-iron supply dynamics [70].
Results for carbon biomass density of fish consumers greater than 10 cm (b10) and 30 cm (b30) depict a bottom-up propagation of the physical forcing over time, peaking ∼9-12 months later than large phytoplankton productivity ( figure 3(A)). By construction of b10/b30, where fish species are classified in terms of body size, bottom-up indirect effects of climate anomalies seem to overbalance species-specific direct temperature-recruitment processes [65]. A feasible process for this bottom-up propagation is a recruitment increase of short-lived species caused by enhanced food supply at lower trophic levels [36]. The spatial response of physical-biological variables for lags 0, 9 and 24 is shown in figures 3(B)-(D). At lag 0, cold SST anomalies and enhanced lphy/tsb are present over the Niño3 region ( figure 3(B)). At lag 9, a mature La Niña is present, together with increased b10 over the Niño3 area and lphy/tsb to the north and south. This northward-southward advection of lphy appears related to Ekman transport induced by enhanced trade winds (supplementary material figure S7). At lag 24 increased b10 also appears displaced to the north and south, while b30 is only located to the north.
A broader picture of this evolution is given in figure 4 in the form of Hovmöller diagrams. Enhanced lphy/tsb/tcb tend to overlap with cold SST anomalies, preserving the timing from figure 3(A). Growth starts at small negative lags over the central tropical Pacific, then anomalies propagate to the east until lag 3 (figure 4(A)) and split thereafter into two branches to the north and south ( figure 4(B)). The long persistence of lphy/tsb/tcb anomalies may be caused by enduring presence of nutrients (nitrate and dissolved iron). As opposed to SSTs, surface nutrients do not interact with the atmosphere [8,36]. Lastly, whereas b10 show a similar response to lphy ∼9 months later, increased b30 are constrained over the northeast tropical Pacific from lags 14 to 24 (figures 3(D) and 4(B)). The latter may be caused by regional currents disposition leading to latitudinal asymmetries in large fish abundance (cf contours in figure 1(C) as a hint). Nevertheless, a species-specific analysis (beyond the scope of this letter) is needed to answer this question.
The evolution of physical-biological variables shown in figures 3 and 4 are in line with the shape and progression of total catch and total system biomass predictable regions in figure 2 (large crosses) and supplementary material figures S3 and S4. This is consistent with the propagation in time of climate anomalies within the food web.
To evaluate the role of diazotrophs dynamics in the physical-biological mechanism, the analysis in figure 4 is replicated in the EcoG-fiswodiaz simulation (identical to EcoG-fis but without diazotrophs dynamics, cf table 2). The differences obtained between the two simulations are exiguous and therefore not shown. This seems consistent with the fact that primary production linked to nitrogen fixation is much less compared to total primary production, particularly in highly-productive marine areas [42,71,72].
To test the sensitivity of ecosystem response to top-down fishing control and model setup, in supplementary material figures S8-S10 we perform the same analyses on EcoOcean and BOATS fished and unfished model output scenarios (  while fishing visibly affects large predators over largescale areas. The differences observed between models (less persistent and westward displaced anomalies of all-size fish biomasses in BOATS model outputs) may be caused by different model complexity, internal hypothesis regarding species dynamics and species relationships, ecosystem functioning and fisheries schemes [42,59]. Lastly, a similar analysis is conducted for Macroecological, although it only provides annual model outputs, does not include a fishing scenario, and propagation in time of anomalies through the food chain is unresolved (supplementary material figure  S11). As expected, the model returns the same picture as in figures 3(B) and (C) for lagged years 0 and 1, with the difference that consumer response (tcb/b10/b30) to the forcing is no longer lagged in time.

Summary and discussion
In this work, the predictability of seasonal to multiannual tropical Pacific marine ecosystem resources was investigated on observations (SAUP catch data) and ecosystem model simulations (FishMIP). Evidence indicates that tropical Pacific fisheries can be forecasted over large-scale areas, and up to three years in advance, by solely considering boreal summer equatorial Atlantic SST anomalies (Atlantic Niños/Niñas) as predictor.
The outstanding skill for predicting Pacific fisheries provided by Atlantic Niños/Niñas can be attributed to the dynamical processes they trigger in the tropical Pacific. These processes involve the alteration of upwelling conditions by anomalous surface winds through the propagation of equatorial Kelvin waves [31,32]. Contrastingly, other well-known remote SST ENSO precursors are associated with more advective processes over the western and central Pacific [19,30] (e.g. NTA, Indian Ocean Dipole), which do not seem to have a substantial impact on highly-productive upwelling-dependent marine ecosystems located further to the east in the basin. These results highlight the importance of equatorial Atlantic impact on ENSO [16,28,68,73].
This study is also supported by previous works that laid foundations for skillful seasonal to multiyear predictability of marine primary productivity in the tropical Pacific [8,36,38]. By providing an extension for these analyses across the entire marine food web and identifying its predictability sources, the results presented here may contribute to a more sustainable management of fisheries in the region.
It remains uncertain whether the described biophysical Atlantic-Pacific relation will persist in the future. It is envisaged that under greenhouse warming conditions biomass over the tropical oceans will decline [59], strong eastern ENSO events will become more likely [74] and the Atlantic-Pacific teleconnection (nonstationary in nature [16,32]) will weaken [58]. Furthermore, future scenarios for fisheries are currently even more uncertain than pre-COVID-19 ones [75].
Although disentangling the interplay of all these factors may require a suite of sensitivity experiments beyond the scope of this study, we provide exploratory results for the period 2021-2054 in supplementary material figure S12 (RCP8.5 fished scenario; cf table 2). Overall, results analogous to those obtained for the historical period are found, supporting the applicability of our findings for the decades to come.
Nevertheless, we simply propound a few possible alternative approaches to marine ecosystem modeling. When additional FishMIP simulations become available, a more comprehensible picture may unfold.

Data availability statement
Annual catch data from FAO major fishing areas are available from the Sea Around Us Project (www.seaaroundus.org/data/#/fao). Hadley Centre Sea Ice and Sea Surface Temperature data set (HadISST) are available from the Met Office Hadley Centre observations datasets (www.metoffice.gov.uk/hadobs/hadisst/). National Centers for Environmental Prediction (NCEP) reanalysis data are available from the NOAA/OAR/ESRL PSL website (https:// psl.noaa.gov/). GFDL-COBALT and FishMIP simulation data (EcoOcean, BOATS and Macroecological) are accessible through the Potsdam-Institute for Climate Impact Research Earth System Grid Federation (ESGF) data node (https://esg.pik-potsdam.de/ search/isimip/).

Acknowledgments
This research was funded by the EU H2020 project TRIATLAS (No. 817578), the Universidad Complutense de Madrid project FEI-EU-19-09 and the Spanish Ministry of Economy and Competitiveness project PRE4CAST (CGL2017-86415-R). This work also acknowledges the 'Severo Ochoa Centre of Excellence' accreditation by the Spanish Ministry of Science and Innovation (CEX2019-000928-S) to the Institute of Marine Science (ICM-CSIC). We thank Derek Tittensor (UNEP-WCMC) and Iliusi Vega (PIK-Postdam) for help with FishMIP data extraction. We thank Roberto Suárez-Moreno (LDEO-Columbia University), Jeroen Steenbeek (ICM-CSIC) and Charles Stock (NOAA-GFDL) for their cooperation and help during the progress of this study. Finally, we would like to thank the two anonymous reviewers for their helpful comments and suggestions, which contributed to improve this manuscript.

Code availability
S4CAST originally published MATLAB® code is open-access and available from the Zenodo repository (doi:10.5281/zenodo.15985) in the URL https:/ /zenodo.org/record/15985. The rest of MATLAB scripts used in this analysis can be made available upon request from the corresponding author (i.gomara@ucm.es). Updated S4CAST versions can be requested to the TROPA-UCM research group (brfonsec@ucm.es).