European small pelagic fish distribution under global change scenarios

The spectre of increasing impacts on exploited fish stocks in consequence of warmer climate conditions has become a major concern over the last decades. It is now im-perative to improve the way we project the effects of future climate warming on fisheries. While estimating future climate-induced changes in fish distribution is an important contribution to sustainable resource management, the impacts on European small pelagic fish—representing over 50% of the landings in the Mediterranean and Black Sea between 2000 and 2013—are yet largely understudied. Here, we inves-tigated potential changes in the spatial distribution of seven of the most harvested small pelagic fish species in Europe under several climate change scenarios over the 21st century. For each species, we considered eight Species Distribution Models (SDMs), five General Circulation Models (GCMs) and three emission scenarios (the IPCC Representative Concentration Pathways; RCPs). Under all scenarios, our results revealed that the environmental suitability for most of the seven species may strongly decrease in the Mediterranean and western North Sea while increasing in the Black and Baltic Seas. This potential northward range expansion of species is supported by a strong convergence among projections and a low variability between RCPs. Under the most pessimistic scenario (RCP8.5), climate-related local extinctions were expected in the south-eastern Mediterranean basin. Our results highlight that a multi-SDM, multi-GCM, multi-RCP approach is needed to produce more robust ecological scenarios of changes in exploited fish stocks in order to better anticipate the economic and social consequences of global climate change.


| INTRODUC TI ON
European Small Pelagic Fish (SPF) species have a key economic and ecological role Fréon et al., 2005). SPFs are the main harvested fish group worldwide, representing between 20% and 30% of the global commercial landings, depending on yearly environmental fluctuations and fishing effort (FAO, 2011). These species accounted for 17% of the E.U. catches in 2015 (European Commission & DG-MARE, 2018) and up to 53% of the Mediterranean and Black Sea landings between 2000 and 2013, but most of them are significantly impacted by overfishing (FAO, 2016).
They are also the main prey of most piscivorous fishes, cetacean and seabirds (Bachiller & Irigoien, 2015;Cury et al., 2011), transferring organic matter from the base of the food web to upper trophic levels (Cury, 2000). Abundantly found in upwelling ecosystems (Cury, 2000), SPFs cover a wide range of regions, therefore overlapping with a large range of environmental conditions ). European Seas, especially the Mediterranean Sea host a high diversity of SPFs, including both temperate-cold and temperatewarm water species Coll et al., 2010).
In the context of climate change, a potential range shift of SPFs may induce (a) major economic and social consequences-especially for countries that rely on fisheries for protein supply (Tacon & Metian, 2009)-and (b) deep changes in ecosystem and food web functioning (e.g. Chaalali et al., 2016). Predicting climate-induced range shifts of these largely harvested species is therefore essential for sustainable resource management and food security (Cheung et al., 2013).
Contemporary fisheries management is mainly based on stock assessment models (e.g., Methot & Wetzel, 2013) that estimate populational parameters (e.g., recruitment, abundance and sustainable harvesting rate) for a given fish stock. Over the last decade, several modelling procedures were developed to include environmental variability in fish stock assessment by focusing on the environment-recruitment relationship (e.g., MacKenzie et al., 2008;Tommasi et al., 2017). Considering environmental uncertainty in stock recruitment is a preliminary step towards a possible integration of global climate change impacts on stock dynamics, however (Edgar et al., 2019;Lee et al., 2018;Punt et al., 2014). In the context of ecosystem-based fisheries management (EBFM), several modelling techniques have been developed to integrate trophic, environmental and societal factors (e.g., Romagnoni et al., 2015;Smith et al., 2015) in order to thoroughly evaluate the status of fish stocks (Forrest et al., 2015). One caveat is that robust stock assessment or EBFM requires data-intensive models (Hilborn, 2011;Pauly, 2000;Wetzel & Punt, 2011), an approach not yet applicable to routine use at large-scale and long-term (Edgar et al., 2019). There is therefore an urgent need for basin-scale longterm management plans that include the combined effects of fishing pressure and climate change, for effective conservation of SPFs and sustainable fisheries management (Faillettaz et al., 2019).
Over the last decades, species distribution modelling has been intensively used to project the effects of past and future climate change on the distribution of suitable habitat for species of conservation concern. (Beaugrand et al., 2019; Bellard et al., 2013;Cheung et al., 2009). Species Distribution Models (SDMs), statistical tools based on the niche-biotope duality (Colwell & Rangel, 2009;sensu Hutchinson, 1978) to conceptualize and investigate biogeographical patterns in relation to environmental conditions, are a popular way to assess which species will be under most threat in a near future and/or which regions will be the most impacted by a reorganization of communities (Sinclair et al., 2010). While stock assessments are based on parameters related to stock dynamics (e.g. spawning biomass, recruitment, growth, mortality), these techniques model the ecological niches of species using (a set of) environmental conditions where the species has been observed. The modelled niche can then be used to project potential distribution of species under different environmental conditions on a broad temporal and spatial scale (Colwell & Rangel, 2009), allowing the investigation of potential future range shifts. The capacity of SDMs to produce long-term, large-scale and comparable (i.e. between species, regions or management zones and periods) future projections is of major importance for species conservation and management (Hollowed et al., 2013). The main objective of SDM-based approaches is the production of robust scenarios of future species distribution for reliable management and conservation perspectives (Cheung, Frölicher, et al., 2016;Goberville et al., 2015;Stock et al., 2011) and best practices recommend multimodel ensemble projections (Buisson et al., 2010), that is the use of a large set of SDMs and climate models  Wilby & Dessai, 2010). Most SDM applications on European SPFs did not include long-term distributional range projections at the regional or European scales (Brown et al., 2006;Maynou et al., 2020;Sabatés et al., 2006;Tsikliras, 2008), however, and a few studies addressed this challenge at the European scale Raybaud et al., 2017) but only for a few SPFs.
Our study aims to address this gap in knowledge by examining (a)  , 1982eneral Assembly, 1982 and discussed possible economic consequences of climate change on the allocation of fishing effort.

| Description of the modelling framework
We used the future Environmental Suitability Index (ESI; i.e. spatialized index ranging from 0 to 1, based on suitability estimated from contemporary 1990-2017 conditions) of the seven SPFs, retrieved from Schickele et al. (2020). This recently developed modelling framework includes (a) a spatial and environmental sampling bias reduction, (b) the use of the convex hull method to generate pseudo-absence, (c) a numerical and ecological evaluation of model outputs and (d) the quantification of uncertainties associated to the selection of SDMs. Environmental parameters used to assess contemporary species distribution (Table 1), the calibration procedure and how we selected the most accurate models (Table 2) are therefore only briefly discussed in this section.
Contemporary  distributions of SPFs were obtained using an ensemble forecasting framework that select-among eight different statistical algorithms-the models that best reproduce the observed spatial distribution of each species (Araújo & New, 2007;Buisson et al., 2010;Pearson et al., 2006). To account for the source of uncertainty related to the choice of a given species distribution model, we considered seven algorithms computed from the Biomod2 package (Thuiller et al., , 2016 itic areas (i.e., >50 km from the coast, independently of depth) or in low-bathymetric (i.e., between 0 and 300 m depth) oceanic regions (Schickele et al., 2020). To prevent from multicollinearity and unnecessary model complexity (Dormann et al., 2007), only one parameter from clusters of correlated parameters (Pearson's r correlation >0.7) was retained (Leroy et al., 2016).
For each species and combination of environmental parameters, the calibration data set was then filtered in an environmental space to reduce sampling bias as much as possible (Varela et al., 2014).
Pseudo-absences were then generated in the same filtered environmental space outside the corresponding convex hull of observation (i.e. considered as a proxy of environmental suitable conditions; Cornwell et al., 2004;Getz & Wilmers, 2006), excluding the 2.5 and 97.5 percentiles. Finally, we used the Continuous Boyce Index (CBI; Hirzel et al., 2006)-the appropriate evaluation metric for a presence/ pseudo-absence calibration data set-to evaluate the robustness of model outputs (see discussion in Leroy et al., 2018): each model with a CBI value over 0.5 was retained (e.g., Faillettaz et al., 2019). Finally, we calculated the response curve of each environmental parameters by keeping other parameters at their mean values for modelled species, using the evaluation strip method (Elith et al., 2005). The ecological realism and relevance of the models were then corroborated by an expert-based inspection/validation of each response curve in order to discard spurious responses to environmental factors (e.g. bimodal response to temperature).

Models (GCMs)
To project the future ESI of each species, we considered five GCMs in the literature, that is from 2.6 to 8.5 W/m 2 , we used three RCPs scenarios: (a) the optimistic peak and decline (RCP2.6), the intermediate "stabilization" (RCP4.5) and the "business as usual" (RCP8.5) scenarios . We considered SSS as constant over time because its temporal variance is negligible in comparison with its spatial variance (Dickson et al., 1988;Faillettaz et al., 2019;Le Marchand et al., 2020). While the spatial variance of SSS allowed us to discriminate marine from brackish waters (e.g., from 35 to 15 between the west and the east of the Danish strait), its temporal var-

| Pre-treatment of future temperature data
Because temperature-related parameters (Table 1) (Table 1). For each GCMs, RCPs, and geographical cell, we therefore estimated the difference between the two datasets and corrected the model-based temperature data accordingly. This process, already applied by Péron et al. (2012) and Cristofari et al. (2018), ensured (a) a perfect correlation (Pearson coefficient r = 1), (b) no RMSD and (c) the same SD between the two data sets for a common period. Results from the correction procedure and corresponding anomalies are shown in Supplementary Material 3.

| Projection of future environmental suitability
Projections of ESI values were carried out at spatial resolutions suitable for either ecological analyses, that is on a 0.1° × 0.1° spatial grid obtained from a linear spatial interpolation (Goberville et al., 2015), or Bogue policy application, that is in each Exclusive Economic Zone (Zeller et al., 2016). Among the existing fishery zones, we chose to focus on EEZsarea stretching from the coastline out to 200 nautical miles over which a country has special rights regarding the use of marine resources-as they are (a) basic units for fisheries management (e.g. attribution of maximum allowed catches by EEZs) and conservation perspectives (Allison et al., 2009;Zeller et al., 2016), (b) at a spatial resolution well-adapted for biogeographic research (Claus et al., 2014) and (c) commonly used in the literature to project the socioeconomics consequences of climate change on fisheries Sumaila et al., 2011Sumaila et al., , 2015

| Uncertainties related to future projections
The selection of a GCM may greatly influence the projected distributions of a species (Goberville et al., 2015): GCMs may diverge for technical or parameterization reasons, may simulate ocean-atmosphere processes in different ways or may vary due to their initial spatial resolution (Beaumont et al., 2008;Goberville et al., 2015;Wiens et al., 2009).
Because of their wide variety and complexity, and because we cannot identify a model that performs better than another (Martinez-Meyer, 2005), it is essential to consider a full range of GCMs to examine the full range of potential future species distributions (e.g. Friedlingstein et al., 2013;Shepherd, 2014). In our study, we adapted the ensemble modelling for future projections: for each set of environmental parameters and each statistical algorithm, five GCMs and three RCP scenarios were considered to project future environmental suitability. For a given period and for each RCP, we performed 10 cross-validation runs, leading to the production of 50 simulations (5 GCM × 10 cross-validation runs) per statistical algorithm. We then computed the corresponding SD among SDMs (i.e. the variability related to the calculation of the ecological niche) and GCMs (i.e. the intrinsic variability linked to the climate system and expected climate conditions) to fully explore the uncertainty related to future environmental suitability projections (Goberville et al., 2015).

| Future environmental suitability
Here, we present for each of the seven SPFs, the projected ESIs in the spatial domain ranging from 10 to 70°N and −30 to 45°E. Species distributional range under RCP8.5 conditions for the late 21st century (2090-2099) are detailed in Figure 1 while other RCPs and periods are provided in Supplementary Material 5.

| Temperate-cold water species
Strong northward shifts in the distribution of ESI are expected for "temperate-cold" water species (Figure 1;

| Temperate-warm water species
Expected changes in ESI values for "temperate-warm" water species (Figure 1;

| Uncertainties in expected environmental suitability
For all SPFs, our projections showed only weak variations, with mean SD ranging from 0.1 to 0.4 (Figure 1). This demonstrates the spatial convergence of our simulations based on a multi-SDM and multi-

| Climatic range shifts between Exclusive Economic Zones
Here, we explored the consequences of potential distribution shifts at the scale of EEZs (Figure 2), manageable units commonly used for projecting the possible socioeconomic impacts of climate change on fish stocks . For each EEZ and SPF,

| Temperate-cold water species
The four temperate-cold water species are of major importance in European fisheries (FAO, 2020), especially along the Atlantic façade (c.a. 50,000 t/year per EEZ), except for European sprat that is largely harvested in the Baltic Sea (Figure 2). These species are currently mostly captured in regions where we found high ESI values over the period 1990-2017 (Figure 2; Schickele et al., 2020). Our models projected a decrease in the ESI values in southern and western Europe (Figure 2). While the Mediterranean EEZs showed a decrease in mean ESI values from 0. 48 (1990-2017) to 0.39 (RCP2.6) or 0.24 (RCP8.5), we projected an increase in the EEZs of the Baltic Sea from 0. 43 (1990-2017) to 0.51 (RCP2.6) or 0.59 (RCP8.5; Figure 2). Our results highlight that a potential mismatch between current fisheries areas and changes in the species climatic range of temperatecold water species could occur by the end of the century; a major decrease in ESI values was for example expected in Morocco and Turkey (Figure 2), that is where species are currently abundantly captured. In contrast, we projected that ESI values may remain steadily constant over the current century in Denmark, the main fishing area for European sprat (Figure 2). At the European scale, the absolute variation of ESI values (i.e. relative to 1990-2017) was expected to range from 22% (RCP2.6) to 33% (RCP8.5), suggesting a potential reallocation of temperate-cold water species population in fisheries management zones.
The projected decrease in ESI value in the Mediterranean Sea, as a response of the expected increase in temperature, may be explained by a limitation of eco-physiological processes (Dahlke et al., 2020) that hinders the survival and/or development of SPFs (Lehodey et al., 2006;Perry et al., 2005;Torri et al., 2018). effects on egg quality (Dahlke et al., 2020). While sea water salinity may influence the specific gravity of marine pelagic fish eggs, that is their vertical distribution in the water column, egg size of SPFsthat has a fundamental impact on the capacity of young larvae to be active, grow, and survive-is under direct influence of SST (Huret et al., 2016;Peck et al., 2013;Torri et al., 2018). SPFs are planktonic feeders (Bachiller & Irigoien, 2015) which depend on high productive areas such as the gulf of Gabès or the eastern North Sea during early development (Rijnsdorp et al., 2009). Despite the projected warming in the Mediterranean Sea, temporal changes in primary production may remain negligible in comparison with its spatial variance: using a multimodel approach, Macias et al. (2015) report a slight increase in primary production of about 9.5 mmol N/m 2 over the period 2015-2095. This supports why primary production did not greatly influence our projections. At an ecosystem scale, SPFs have a pivotal role (Chaalali et al., 2016) in marine food webs, contributing to carbon fluxes from lower trophic level to top predators (Cury, 2000).
Projected distributional range shifts-and potential ensuing changes in abundance patterns (Helaouet & Beaugrand, 2009;Kulhanek et al., 2011;VanDerWal et al., 2009)-may deeply modify the Mediterranean or Baltic seas ecosystems through trophic cascading effects on the upper trophic levels which feed on these species, including negative effects on their fisheries (Maynou et al., 2014). This climatic resilience-related issue may be assessed through the development and adaptation of ecosystem management and protection strategies (Giakoumi et al., 2017;e.g. McLeod et al., 2009), in regions or for species (e.g. SPFs) identified as sensitive to climate-induced changes (e.g. the Mediterranean Sea; Figure 2).

| Management and economic implications
While fishing has expanded into the high seas over the last decades as a result of an increasing demand for fish and the overexploitation of coastal waters (Sumaila et al., 2015), quantifying changes in the allocation of fisheries catches by maritime country and within EEZ waters allows to focus on a spatial scale that is politically and economically viable Zeller et al., 2016).
Depending on the EEZ, SPFs may experience high fishing mortalities, especially in the Mediterranean Sea (i.e. twice the maximum sustainable yield; FAO, 2016), with putative negative impacts on species growth, reproduction and stock production (Brander, 2007;Fréon et al., 2005;Lehodey et al., 2006). In our study area, 65% of the

| Perspectives on modelling small pelagic fisheries
Despite their popularity, SDMs have inherent limitations-depending on their application context-such as the assumption of niche conservatism (Peterson & Soberón, 2012) and the failure to integrate species interactions or dispersal processes (Araújo & Guisan, 2006). While recent advances in species distribution modelling proposed to overcome such shortcomings (e.g. dispersal constrained SDM; Boulangeat et al., 2012;e.g. joint-SDMs;Harris, 2015), most of these ecological processes are mainly associated with local changes (Beaugrand & Kirby, 2018). These perspectives are at the cost of the amount of data needed to calibrate the models which can impede their broad scale-applicability, one of the main strength of SDM (Marmion et al., 2009). Based on our multi-GCM, multi-SDM and multi-RCP approach that integrated climate uncertainty (Friedlingstein et al., 2013;Goberville et al., 2015;Shepherd, 2014) and to contribute to increase the probability of success in fishery management strategies (Jones et al., 2012), we encourage further research in the regions that we have identified as the most vulnerable, such as the Mediterranean Sea. The whole Mediterranean Sea ecosystem, including highly impacted regions (Bănaru et al., 2013;Hattab et al., 2013;Piroddi et al., 2015), has been extensively studied over the last decades, and it is now well documented that changes in small pelagic populations-due to fishing or natural drivers-will strongly alter ecosystem structure and functioning (Palomera et al., 2007). Considering results from SDMs in combination with ecosystem models for fisheries management will help to better anticipate the consequences of climate-induced distributional shifts in small pelagic fish on the whole ecosystem (Chaalali et al., 2016). The need is urgent as many countries in the Mediterranean Sea are directly or indirectly dependant on activities that involve exploitation of marine fish resources (Selig et al., 2019).

ACK N OWLED G EM ENTS
This paper is dedicated to the memory of Pr. Patrice Francour who passed away on 13 October 2019. He dedicated his life to the protection and understanding of Mediterranean ecosystems, initiated and co-supervised this work. We acknowledge the Ocean Biogeographic Information System (OBIS), the Global Biodiversity Information Facility (GBIF) and FishBase for providing species occurrence data.
Finally, we thank the IPCC Coupled Model Intercomparison Project (CMIP) and the climate modelling groups for making available their model output.

CO N FLI C T O F I NTE R E S T
The authors have no conflict of interest to declare.

AUTH O R CO NTR I B UTI O N
VR and PF conceived and supervised the study. AS, VR and EG collected the data. BL provided the initial modelling framework. AS performed the numerical analysis. GB, BL, TH, EG and VR helped in the modelling process. AS wrote the first draft. All authors made substantial contributions to the successive versions of the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
These data were derived from the following resources avail-