Climate change impacts on living marine resources in the Eastern Tropical Pacific

Project shifts in the habitat suitability of 505 fish and invertebrate species in the Eastern Tropical Pacific that are likely to occur by the mid‐21st century under “high greenhouse gas emissions” (RCP 8.5) and “strong mitigation” (RCP 2.6) scenarios.

Rapid changes in the ocean are exposing ectotherms to unfavourable conditions Pörtner, 2001) that lower their growth, reproduction and survival rates. These physiological impacts lead to shifts in species biomass and distributions (Cheung et al, 2010;Gattuso et al., 2015;Pörtner et al., 2014). The vulnerability of a species to climate change mainly depends on the scope between current conditions and the species physiological tolerance limits (Pörtner & Peck, 2010).
Tropical species generally exhibit narrow thermal tolerances relative to temperate species and therefore tend to be more sensitive to rapid changes in environmental temperature (Jones & Cheung, 2015;Madeira et al., 2012;Nguyen et al., 2011;Pörtner & Peck, 2010). Under the high-emissions scenario (Representative Concentration Pathway 8.5), species richness in the tropics is projected to decrease by more than 20% by 2050 relative to the 2000s (Jones & Cheung, 2015). Maximum catch potential (a proxy for maximum sustainable yield) is also projected to decrease globally by 3.4 million tonnes per degree Celsius of atmospheric warming, with the tropics being most at risk of such impacts (Cheung et al., 2010. These climate-driven declines in the availability of tropical living marine resources have the potential to affect national economies and food security, especially in developing countries with low adaptive capacity (Allison et al., 2009;Blasiak et al., 2017;Lam et al., 2016;Lam et al., 2020).
Despite the profound impacts that climate change will have on tropical marine ecosystems, most of the available information is for northern, temperate regions (Poloczanska et al., 2016). Given the widespread lack of access to scientific vessels in tropical countries, fishing vessels represent useful platforms to gain information on species distributions (Wehrtmann et al., 2012). Consequently, patterns in species caught by fisheries over time can help us understand regional biogeographic shifts occurring in response to climate change (Maestri, 2020). Within the Eastern Tropical Pacific Ocean (ETP) (Figure 1), climate change effects on fisheries have only been documented for the Gulf of California (Páez-Osuna et al., 2016). However, climate impacts on fishery species may go undetected due to the confounding effects of climate change, climate variability and overfishing (Hsieh et al., 2006), which are hard to disentangle due to the absence of long-term fisheries and oceanographic time series (Pauly & Zeller, 2016). Moreover, while global projections show that, in general, living marine resources will shift polewards (Jones & Cheung, 2015;Poloczanska et al., 2013;Pörtner et al., 2014), observations have underscored the variability in climate change responses across species and space (Lenoir & Svenning, 2015;Morley et al., 2018;Pinsky et al., 2013), with consequences for fisheries management and conservation actions (Arafeh-Dalmau et al. 2020;Frazão Santos et al. 2020).
This study aims to examine the subregional range shift patterns of key fishery species in the ETP and the implications of such shifts for key fisheries and regional biogeography. We hypothesize that climate change will force species to shift towards cooler waters brought up to the surface through upwelling. We use an ensemble of species distribution models to project the future geographic distribution of living marine resources across the ETP under two climate change scenarios. This approach is based on the concept of environmental niche, defined as a set of environmental conditions that allow a species to persist (Hutchinson, 1959).
Species distribution models can help elucidate the complexity of biogeographic responses of marine species at subregional scales (Guisan & Thuiller, 2005;Guisan & Zimmermann, 2000) even in data-limited areas such as the ETP (Hernandez et al., 2006;Wisz et al., 2008). Based on the model projections, we discuss the potential impacts of climate change on resource availability and regional marine biogeography.

| Oceanographic setting of the Eastern Tropical Pacific
The Eastern Tropical Pacific Ocean is defined here as the area between 31 °N and 5 °S, from the northern Gulf of California to northern Peru, with the East Pacific Barrier to the west and the F I G U R E 1 Exclusive economic zones in the Eastern Tropical Pacific of interest in this study. 1: Mexico; 2: Guatemala; 3: El Salvador; 4: Nicaragua; 5: Costa Rica; 6: Panama; 7: Colombia; 8: Ecuador; 9: Galapagos; 10: Peru Central American Isthmus to the east (Figure 1). Upwelling systems within the ETP play an important and complex role in driving the spatial patterns of temperature, primary productivity, pH and oxygen (Fiedler & Lavín, 2017). As such, the ETP is delineated by the California Current eastern boundary upwelling system in the north and the Humboldt Current eastern boundary upwelling system, South Equatorial Current and equatorial upwelling system to the south ( Figure 2; Fiedler & Lavín, 2017 produced by wind jets across Central America (Fiedler & Lavín, 2017). Oxygen minimum zones (OMZs) are formed below the shallow thermocline (<80 m) as a result of the high primary productivity in upwelling systems, strong stratification and sluggish circulation ( Figure 3). The southern ETP is characterized by low salinity (<34 psu), while the northern limits are characterized by high salinity (>34 psu), except for an area off Baja California with low-salinity waters transported by the California Current (Fiedler & Lavín, 2017).  Pauly & Zeller, 2015). Marine fisheries in the ETP are diverse in terms of species caught, gears used and fleet sizes (Lluch-Cota et al., 2018).
This study focused on species caught by four key fisheries in the ETP: small-scale fisheries, shrimp trawl fisheries, small pelagic fisheries and large pelagic fisheries. Small-scale fisheries generally operate in coastal areas within the continental platform, using a wide variety of gears to capture a large diversity of species. In contrast, the fisheries for shrimp, as well as small and large pelagics, are predominantly large scale. Shrimp trawl fisheries mainly target penaeid, pandalid and solenocerid shrimp species. The small pelagic fisheries mainly operate within the continental shelf and target engraulids and clupeids. Large pelagic fisheries generally operate in more offshore oceanic waters and target tuna, mahi mahi, swordfish and sharks (Cisneros-Montemayor et al., 2013;Donadi et al., 2015;Haas et al., 2015;Harper et al., 2014;Lindop et al., 2015;Lluch-Cota et al., 2018;Schiller et al., 2014;Trujillo et al., 2015).

| Biotic data
We modelled species caught by the four major fisheries in the ETP, irrespective of the species commercial value. While all fisheries within the region catch a wide variety of species, only a subset of them is landed and sold (Appendix S1: Suppl. 1). Landings records have a low degree of taxonomic resolution and do not include discards. To obtain an accurate representation of species caught by fisheries in the region, inclusive of bycatch and discards, we conducted a thorough scientific literature review that yielded a list of 652 species, comprised of 512 bony fish, 74 elasmobranch species, 47 crustacean species, 16 mollusc species and three echinoderm species. In this study, we attempted to model all species caught by these fisheries, as low-valued species may become valuable in the future.

| Abiotic data
We applied generalized linear models (GLM-identity link, Gaussian distribution) to statistically bias-correct (Wilby et al., 1998) average annual climatologies for bottom and surface temperature, oxygen, salinity, pH, primary productivity and mixed layer depth for the ETP (Appendix S1: Suppl. 4). Earth system models (ESMs) account for processes in coarse-grid cells that lead to biases at the local scale.
Here, we attempt to correct these biases by modelling the relationship between ESM output and observed climatologies.  (Giorgetta et al. 2013). The performance of these ESMs has been extensively examined and tested for applications to the marine realm (Kwiatkowski et al. 2017;Laufkötter et al. 2015). All environmental parameters were regridded and interpolated on to a global 0.5° longitude × 0.5 latitude raster using the bilinear interpolation method, before the bias correction (see Lam et al., 2016).
To account for spatial autocorrelation, we also included the interaction between latitude and longitude as an independent variable in all models. Depth was included as an independent parameter for models of bottom environmental conditions. We did not bias-correct pH, because time series for observed pH data do not exist.
We used the models to produce a regional annual climatology for each parameter based on outputs for each ESM. The global climatologies used to train the species distribution models include the regional bias correction. We assumed the statistical relationships between ESM and observed climatology will hold in the future and, therefore, used them to project future observed environmental conditions given a set of ESM projections (annual averages of each environmental parameter for RCP 2.6 and RCP 8.5 from 2001 to 2060).

| Species distribution models
The current and future distributions of the 547 focal species were projected using species distribution models (SDMs). We used a multimodel approach (Jones & Cheung, 2015) to account for the variability across different Earth system models and SDM outputs and increase the accuracy of the projections (Guo et al., 2015;Jones & Cheung, 2015). We applied four SDMs to quantify the environmental niche of each species: Surface Range Envelope (Araújo & Peterson, 2012), maximum entropy method (Maxent) (Elith et al., 2011), generalized boosting model (Elith et al., 2008) and artificial neural networks (Lek & Guégan, 1999). The input data for each SDM were the species occurrence raster (75% of the original presence data were used to train the model) and the global climatology of the environmental conditions. We selected variables representing bottom water conditions for benthic and demersal species, and surface water conditions for pelagic or coastal species (Froese & Pauly, 2018;Robertson & Allen, 2015). We avoided overparameterization by selecting the subset of the environmental parameters that resulted in the highest specialization (narrowness of the niche) and marginality (difference between the niche and the available environment) values (Appendix S1: Suppl. 5) produced by the ecological-niche factor analysis (ENFA, Basille et al., 2008;Calenge, 2006).
For each species, we ran the four SDMs using the three global climatologies (GFDL-ESM-2G, IPSL-CM5-MR and MPI-ESM-MR), resulting in outputs from a total of 12 models per species. All SDMs were run with the Biomod2 package in R (Thuiller et al., 2016). Each SDM calculated a Habitat Suitability Index (HSI) value for each spatial cell in the ETP region, ranging between 0 (not suitable) and 1 (very suitable).
We evaluated the accuracy of each model using and area under the curve (AUC) analysis of the receiver operating characteristic (ROC). We used the ROC to compare the fitted HSI with the species occurrence raster reserved for testing the model fit (25% of the original presences). Models with an AUC below 0.5 were eliminated, as predictions were worse than random (Sing et al., 2005).
This analysis was performed with the pROC package in R (Robin et al., 2011). To be consistent with the approximate time frame represented by the occurrence records and for the development of the models, we used the average SDM predictions for 1970 to

2000.
We then projected changes in the geographic distribution of the species environmental niche for each year between 2001 and 2060 under the "high emission" (RCP 8.5) and "strong mitigation" (RCP 2.6) scenarios. We built an ensemble of model outputs for each species and climate change scenario. Specifically, for each ESM, we first calculated the average HSI weighted by the AUC values of each species distribution model and then averaged HSI values across ESMs, to produce one HSI value per cell. If the habitat suitability was higher than the species prevalence (i.e. the fraction of cells in which the species was present), we considered the species to be present in the cell (Phillips et al., 2009). We averaged projections over a 20-year time frame to reduce the effect of interannual variability of climatic conditions on species distributions (Stock et al., 2011).
We calculated indicators of biogeographic shifts for the species assemblages caught by each of the four fisheries. We present the results for the shrimp trawl fishery separated in target and bycatch species. These indicators include the shift in geographic and depth centroids, local species loss rates, local invasion rates, species turnover (Cheung et al., 2009)  Latitudinal and longitudinal shifts were estimated as the shortest distances between the centroids according to the Haversine method (Hijmans et al., 2018), which assumes a spherical earth. Local invasion and local loss rates were estimated as where n is the number of species per cell at the beginning of the century, and n I i,y and n E i,y are the number of species invading or going extinct in each cell, respectively, by the end of the study's time frame. Finally, species turnover is the sum of invasion and extinction rates.
By the mid-21st century, surface temperatures are projected to increase by 0.96°C, on average, relative to the present, under RCP 8.5, except in upwelling regions. Surface oxygen follows the same spatial patterns, which is not surprising considering the negative correlation between temperature and oxygen solubility.
Amongst the 547 species with sufficient data for species distribution modelling, 505 species were used for further analysis because at least three of the four species distributions models had AUCs above 0.5 (Appendix S1: Suppl. 6). Variability in AUC values was high across SDMs and low across Earth system models (Appendix S1: Suppl. 7

| Species shifts
In the northern region (>15° N), the centroids of species distributions were projected to shift towards the north-west at an average rate of 71 km decade -1 , while in the northern equatorial region (0° to 15° N) and the southern ETP (0° to 20° S), species were projected to shift south-east at an average rate of 59 and 30 km decade -1 , respectively ( Figure 4). The direction of projected species geographic shifts was similar across fisheries. Demersal species in most countries were found to move towards shallower waters by an average rate of approximately 1-13 m decade -1 ; however, there was considerable variability in the direction and magnitude of depth shifts across species ( Figure 5).

| Habitat suitability
The spatial patterns of changes in habitat suitability by 2041-2060 were similar between RCP 2.6 and RCP 8.5, although in most cases the magnitude of change was greater for RCP 8.5 (Table 1) (Figure 8). Projected turnover rates for species caught by small-scale fisheries were higher than for any other fishery, surpassing 38% between Costa Rica and Galapagos, and reaching almost 80% in Colombia (Table 2).
Projected decreases in habitat suitability for species caught by large pelagic fisheries were highest between Guatemala and Nicaragua (−20% to −26%) (Table 1), and projected turnover rates were highest between Guatemala and Ecuador (Table 2). Most invasions were projected to occur from Panama to Peru (Figure 8), while local losses were expected to remain high throughout the study region, except along the northern and southern limits ( Figure 7).
Models projected large declines in habitat suitability for species caught by small pelagic fisheries from Guatemala to Colombia (Table 1) and a 17% increase for Peru (Table 1). Species turnover projections were below 30% (Table 2), with the highest local losses expected along the continental shelf, especially from Costa Rica to Ecuador (Figure 7). Projected invasion rates were also highest for the continental shelf area (Figure 8).
The habitat suitability of species targeted by the shrimp trawl fishery was projected to increase in all EEZs except Nicaragua, and particularly in Mexico and Peru (Table 1, Figure 6). Local loss rate projections were low throughout the study region apart from small areas in the Gulf of California (Figure 7). Shrimp invasion rates were projected to follow a patchy distribution, mainly along the coastline of northern Mexico and southwards of Panama (Figure 8). Target species in this fishery also showed the lowest turnover rates, with projections below 6% for all EEZs (Table 2).
In contrast to species targeted by the shrimp trawl fishery, the habitat suitability of species caught as bycatch was projected to  (Figures 6 and 7). Invasion rates were projected to be much lower than for small-scale fisheries and were limited to the northern and southern limits of the study area ( Figure 8).
Most species in each EEZ were not projected to undergo declines in habitat suitability (Appendix S1: Suppl. 8

| D ISCUSS I ON
Our results provide insights into the expected climate change impacts on regional biogeography and living marine resources in the ETP between now and 2041-2060, a time frame deemed relevant to inform climate change adaptation of fisheries management actions.
Results show limited divergence in the projected oceanographic conditions and habitat suitability losses between greenhouse gas emissions scenarios from the present day to 2041-2060 because of lagged responses of some oceanic variables to changes in atmospheric greenhouse gas concentrations.
Two general patterns emerge when examining the oceanographic processes responsible for the distribution of fisheries species: species are shifting towards cooler waters in the northern and southern margins of the ETP and towards more oxygenated, shallower waters. Temperature and oxygen were predicted by our models to be F I G U R E 6 Projected change in the habitat suitability of species groups caught in the four main fisheries (large pelagics, small pelagics, shrimp trawl (incl. of bycatch) and small scale) by 2041-2060relative to present (2001-2020. Warm hues in the colour ramp denote losses in habitat suitability, while cool hues denote gains the most important variables shaping the environmental niche of the study species, and therefore, warming and deoxygenation will likely drive the redistribution of species in the ETP. These patterns may be representative of the broader response of marine biodiversity to climate change in the region. Such patterns also agree with expectations from proposed theory explaining the biological responses of marine fishes and invertebrates to changing temperature and oxygen levels. For example, the "oxygen-and capacity-limited thermal tolerance" theory (OCLTT) suggests that temperatures above an organism's thermal tolerance threshold results in a smaller aerobic scope for physiological functions, like growth and reproduction (Pörtner, 2001).
According to the Gill-Oxygen Limitation Theory (GOLT), an increase in oxygen uptake to meet the higher metabolic oxygen demands under ocean warming is limited by the constraints of available area for gaseous exchange in the gill . Thus, fish move to waters with temperatures that resemble those of their original habitats to satisfy organisms' oxygen needs.

| Species shifts
The Humboldt Current and California Current eastern boundary upwelling systems (produced by alongshore winds in Peru and northern Mexico) and the equatorial upwelling systems (near the equator produced by the Coriolis force) have a cooling effect along the northern and southern limits of the ETP (Fiedler & Lavín, 2017).
These upwelling systems are separated by the eastern Pacific warm pool along Central America, resulting in an inverse temperature gradient in the Northern Hemisphere (Fiedler & Lavín, 2017).
Consequently, between 0° N and 15° N species are moving towards the equator instead of the poles (Fiedler & Lavín, 2017;Pörtner et al., 2014). As expected, species are shifting at a faster rate in tropical areas with weaker latitudinal temperature gradients, where their preferred temperature has shifted further away. In contrast, the steeper temperature gradient along the Humboldt Current allows for species to find their preferred temperature within shorter distances (Robinson et al., 2015).
Species within the ETP shift towards the equator when waters warm during El Niño events (Sielfeld et al., 2010), further supporting the projected equatorwards shift of tropical species. Local losses, population recoveries and range extensions are common responses to El Niño-Southern Oscillation-related interannual temperature variability in the ETP (Mora & Robertson, 2005). Reports from as early as 1982 show that tropical shrimp shift southwards towards Peru during El Niño (Barber & Chavez, 1983). A more recent study identified 100 tropical species in Chilean waters (with subtropical and temperate climatic conditions) during El Niño years (Sielfeld et al., 2010).
In contrast to observed range shifts elsewhere (Dulvy et al., 2008;Poloczanska et al., 2013), our findings show that species are projected to shift towards shallower instead of deeper waters. This shoaling of species can be attributed to the expansion of oxygen minimum zones (OMZs), which drive most organisms into shallower and more oxygenated waters (Stramma et al., 2008(Stramma et al., , 2010(Stramma et al., , 2012. These shallower waters are also warmer, which increases the energy required to meet basic metabolic demands and may require organisms to compensate for temperatures outside their tolerance range (Gallo & Levin, 2016;Pörtner et al., 2014). OMZs are known to compress the habitats of both benthic and pelagic species (Gallo & Levin, 2016). For example, OMZs have been shown to compress the habitat of billfish in the eastern Pacific warm pool (Prince and Goodyear 2007) and of small pelagics in Peru (Bertrand et al., 2011). The expansion of OMZ was also found to force echinoderms in the California Current to shoal, while the contraction of OMZ during El Niño temporarily expanded the habitat of demersal fish towards deeper waters in Peru (Arntz et al., 2006;Keller et al., 2015;Sato, Levin, and Schiff, 2017).
Decreases in the habitat suitability and local losses of the living marine resources focused on in this study mainly occurred across the EEZs of Central America and Colombia, coinciding with the TA B L E 2 Projected species turnover (%) by 2041-2060 relative to 2001-2020 for species groups caught in the four main fisheries (large pelagics, small pelagics, shrimp trawl (incl. of bycatch) and small scale) in the Pacific Exclusive Economic Zones from Mexico to Peru warming and expansion of the eastern Pacific warm pool (Fiedler & Lavín, 2017). High invasion rates along the northern limits of the study area and south of 10˚N could be caused by species shifting towards the cooler waters of the upwelling systems. Projections for upwelling systems under climate change, however, are uncertain as Earth system models do not resolve upwelling processes well . The potential increase in upwelling intensity could affect species biogeography. For example, enhanced hypoxia and acidification associated with upwelling activity could limit the beneficial effects of cooling and higher primary productivities on habitat suitability (Bakun, 2017;Bakun et al., 2015;Fiedler & Lavín, 2017).

| Implications for fisheries and conservation
Our findings highlight the importance of local-scale oceanographic and biological data to elucidate the multidimensional biogeographic shifts on fishery species and their potential impacts on fisheries in the region. Overall, changes in the habitat suitability, and therefore composition of species caught by the four main fisheries, are expected to be most severe along Central America, with substantial variations in the magnitude of impacts across fisheries.
The results suggest that shrimp fisheries may benefit from the impacts of climate change in the ETP because of the increase in habitat suitability of its target species. Shrimp in general may be less vulnerable to climate change because of their fast population growth rates, high larval dispersal rates and low ecological specificity (Hsieh et al., 2006). Nevertheless, many shrimp stocks throughout the region are overfished (Cisneros-Montemayor et al., 2013;Donadi et al., 2015;Haas et al., 2015;Harper et al., 2014;Lindop et al., 2015;Schiller et al., 2014;Trujillo et al., 2015).

| Model robustness and uncertainty
While the regional geographic patterns of projected changes are considered reliable, the specific magnitude of projected changes is affected by model and scenario uncertainties Wabnitz et al., 2018). While the multimodel approach accounts for uncertainties associated with future emissions, species distribution models' and Earth system models' structures, several additional sources of uncertainty remain (Peterson & Soberón, 2012;Stock et al., 2011). First, the results depend on the accuracy of the biascorrected ESM projections  and the assumption that the correlation between the ESMs and the observed climatology will hold in the future. In addition, model outputs rely on ESMs, which exclude small-scale processes that would allow to better resolve regional upwelling dynamics and coastal processes. Second, the species distribution modelling approach applied here focuses on habitat suitability in a single-species context. Therefore, the method does not account for ecosystem effects, such as interspecies dynamics, the ability of invading species to establish themselves in new habitats, or the possible biodiversity and habitat loss that invading species may cause (Pecl et al., 2017). Local losses and invasions have the potential to modify ecosystem structure and function (Marzloff et al., 2016;Pecl et al., 2017). Third, we may have underestimated the declines in habitat suitability by modelling the realized niches of species rather than populations, for which there are little data available. The average niche breadth of a species is much wider than that of a population and therefore will be less sensitive to projected changes in environmental conditions. Fourth, several tropical species are already close to their thermal limits, likely increasing their vulnerability to further warming (Sunday, Bates, and Dulvy 2011).
Fifth, we do not consider other human stressors that have and will continue to drive impacts on living marine resources (Galbraith et al., 2017;Peterson & Soberón, 2012;Polidoro et al., 2012). For example, our models do not account for the environmental impacts of shrimp fisheries, nor for the effects of overfishing of large pelagics in the ETP (Dent & Clarke, 2015;Espinoza et al., 2018). Therefore, billfish, tunas and other large pelagics are probably even more vulnerable to climate change than our results indicate. Moreover, local losses or invasions of top predators may trigger trophic cascades that cause disproportionate impacts on marine ecosystems (Heupel et al., 2014). Sixth, the possibility of genetic and phenotypic plasticity that determine potential acclimation and rapid evolution of marine species may reduce their sensitivity to climate change (Calosi et al., 2016;Quintero and Wiens, 2013). Evidence of evolutionary responses of marine fishes and invertebrates to climate change is still limited and is an important area for further exploration in future studies. Overall, most of these sources of uncertainty would likely produce stronger declines in habitat suitability, while incorporating the potential adaptation responses of species may result in more optimistic projections.
Our findings contribute to the understanding of species responses to the growing threat of climate change in a complex oceanographic region, in support of efforts to implement management and conservation actions. This study can help identify which species may need greater protection in the future and identify areas that would support greater resilience in the face of climate change, such as biodiversity refuges (Hoffmann, Irl, and Beierkuhnlein, 2019;Kujala et al., 2013). The identification of areas that may be particularly vulnerable to climate change may also be used to inform marine spatial planning for climate adaptation initiatives in the ETP ). Thus, species projections can inform policy decisions and conservation strategies that ensure the protection and sustainable use of living marine resources (Wilson et al., 2020).

DATA AVA I L A B I L I T Y S TAT E M E N T
We used open access data sources for the analysis in this manuscript and have included them in the references and supplements. Tayler aims to understand how environmental change may transform the biogeography of living marine resources, using empirical data and model simulations. She is also developing a new method that will allow us to understand and project the impacts of warming and deoxygenation on data-deficient fisheries resources.

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found online in the Supporting Information section.