Current and Future Influence of Environmental Factors on Small Pelagic Fish Distributions in the Northwestern Mediterranean Sea

Pennino, Maria Grazia; Coll, Marta; Albo-Puigserver, Marta; Fernández-Corredor, Elena; Steenbeek, Jeroen; Giráldez, Ana; González, María; Esteban, Antonio; Bellido, José M.

doi:10.3389/fmars.2020.00622

ORIGINAL RESEARCH article

Front. Mar. Sci., 24 July 2020

Sec. Marine Fisheries, Aquaculture and Living Resources

Volume 7 - 2020 | https://doi.org/10.3389/fmars.2020.00622

This article is part of the Research Topic Managing for the Future: Understanding the Relative Roles of Climate and Fishing on Structure and Dynamics of Marine Ecosystems View all 22 articles

Current and Future Influence of Environmental Factors on Small Pelagic Fish Distributions in the Northwestern Mediterranean Sea

$\r\nMaria Grazia Pennino,,*&#x;$ Maria Grazia Pennino^1,2,3*†

Marta Coll^2,4,5†

Marta Albo-Puigserver⁴

Elena Fernández-Corredor⁴

Jeroen Steenbeek^4,5

Ana Giráldez⁶

María González⁶

Antonio Esteban⁷

José M. Bellido^2,7

¹Centro Oceanográfico de Vigo, Instituto Español de Oceanografía, Pontevedra, Spain
²Statistical Modeling Ecology Group, Departament d’Estadística i Investigació Operativa, Universitat de València, Valencia, Spain
³Departamento de Ecologia, Fishing Ecology Management and Economics (FEME) – Universidade Federal do Rio Grande do Norte – UFRN, Natal, Brazil
⁴Institut de Ciències del Mar (ICM-CSIC), Barcelona, Spain
⁵Ecopath International Initiative Research Association, Barcelona, Spain
⁶Centro Oceanográfico de Málaga, Instituto Español de Oceanografía, Fuengirola, Spain
⁷Centro Oceanográfico de Murcia, Instituto Español de Oceanografía, San Pedro del Pinatar, Spain

In the Northwestern Mediterranean Sea, the European sardine (Sardina pilchardus) and the European anchovy (Engraulis encrasicolus) are the most important small pelagic fish in terms of biomass and commercial interest. During the last years, these species have experimented changes in their abundance and biomass trends in the Northwestern Mediterranean Sea, in addition to changes in growth, reproduction and body condition. These species are particularly sensitive to environmental fluctuations with possible cascading effects as they play a key role in connecting the lower and upper trophic levels of marine food webs. It is therefore essential to understand the factors that most profoundly affect sardine and anchovy dynamics. This study used a two-step approach to understand how the environment influences the adult stages of these species in the Northwestern Mediterranean Sea. First, we explored the effects of environmental change over time using Random Forests and available datasets of species occurrence, abundance, biomass and landings. We then applied species distribution models to test the impact of the extreme pessimistic and optimistic Intergovernmental Panel on Climate Change (IPCC) pathway scenarios, and to identify possible climate refuges: areas where these species may be able to persist under future environmental change. Findings from the temporal modeling showed mixed effects between environmental variables and for anchovy and sardine datasets. Future pathway projections highlight that both anchovy and sardine will undergo a reduction in their spatial distributions due to future climate conditions. The future climate refuges are the waters around the Rhone River (France) and the Ebro River (Spain) for both species. This study also highlights important knowledge gaps in our understanding of the dynamics of small pelagic fish in the region, which is needed to progress towards an ecosystem approach to fisheries management.

Introduction

Small pelagic fish species (SPF) are sensitive to environmental fluctuations, which can be amplified from climate variability, with cascading effects up and down the food web, due to their high biomass at intermediate levels and their potential wasp-waist control of the food web (Cury and Roy, 1989; Cury et al., 2000; Bakun, 2006; Checkley et al., 2009; Fauchald et al., 2011). In the Mediterranean Sea, SPF species such as European sardine (Sardina pilchardus, Walbaum, 1792; hereafter sardine), and European anchovy (Engraulis encrasicolus, Linnaeus, 1758; hereafter anchovy) have been shown to be key elements in the transfer of energy from lower to higher trophic level organisms (Palomera et al., 2007; Coll et al., 2008; Albo-Puigserver et al., 2016). Historically, they showed an important bulk of biomass and production and played a key ecosystem role due to a strong coupling between the pelagic and the demersal environments (Coll et al., 2008), acting as important prey of pelagic predators such as tunas (Navarro et al., 2017), cetaceans (Gómez-Campos et al., 2011) and pelagic seabirds (Navarro et al., 2009; Cury et al., 2011), and of demersal predators (Recasens et al., 1998; Mellon-Duval et al., 2017; Saraux et al., 2019).

Recent changes in SPF have been observed in the Mediterranean Sea, in parallel with an overall increase of fishing effort, changes in environmental variability and a decline of primary productivity (Piroddi et al., 2017). Specifically, in the Northwestern Mediterranean Sea we have observed declines in biomass and landings of anchovy and, especially important, of sardine in parallel with a spatial expansion of round sardinella (Sabatés et al., 2006; Palomera et al., 2007; Van Beveren et al., 2016). Changes in anchovy and sardine landings and abundances have also been related with increases in fishing impact and recent high rates of exploitation (Palomera et al., 2007; Scientific Technical and Economic Committee for Fisheries [STECF], 2016; FAO, 2018; Ramirez et al., 2018).

These changes have been also linked to environmental shifts that can directly influence annual recruitment, growth and the general condition of organisms, for which several hypotheses have been formulated (e.g., Palomera et al., 2007; Martín et al., 2008; Van Beveren et al., 2014; Brosset et al., 2017; Coll et al., 2019). For example, the impact of higher sea surface temperature (SST) could negatively affect some species, such as sardine that prefers cold waters to reproduce, and positively others, such as round sardinella, that prefers warm waters to reproduce (Sabatés et al., 2006; Palomera et al., 2007; Maynou et al., 2014).

Within this context, and in order to provide an exhaustive insight into the environmental influence on SPF, we developed a two-step approach to (1) explore the effect of the environment from a historical temporal perspective and to (2) spatially predict the distribution of SPF species under different future environmental scenarios to identify possible climate refuges.

The temporal perspective was explored applying Random Forests (RFs) (Breiman, 2001) to different historical datasets regarding presence/absence, abundance, biomass, and landings of (mostly adult) sardine and anchovy, with the aim to test the impact of different environmental parameters and temporal lags.

We based the spatial perspective on the Species Distribution Models (SDM) framework. The various algorithms within this framework link spatial species occurrences to environmental variables with the aim to predict if species are likely to occur in unsampled locations or future time periods (Martínez-Minaya et al., 2018). We used residuals autocovariate Boosted Regression Trees (RAC-BRT) (Elith et al., 2008; Crase et al., 2012; Escalle et al., 2016; Barcala et al., 2020) to test different environmental future pathway scenarios as defined by the Intergovernmental Panel on Climate Change (IPCC, 2013) to identify possible climate refuges. The term climate refuges refers to areas where species retreat to, and can persist, under increasing environmental stress with the potential to re-expand once the stress abates (Bennett and Provan, 2008; Keppel et al., 2012). Such areas can act as sources of re-colonization when environmental conditions improve and often have long-lasting imprints on species distributions. We argue that identifying climate refuges can provide essential knowledge to improve the SPF management by identifying key areas for conservation of these species and potential for re-expansions both today and in the future (Monsarrat et al., 2019).

Materials and Methods

Study Area and Species Data

This study covers Geographical Sub-Areas (GSAs) 06 and 07 of the Northwestern Mediterranean Sea, located, respectively, in Northern Spain (between Cabo de Creus and Cabo de Palos), and the Gulf of Lions (France) (Figure 1). This is a highly productive area within the Mediterranean Sea due to the presence of two large rives (Rhone and Ebro) and a counter-clock oceanographic circulation that brings nutrients to the upper slope and shelf areas (Estrada, 1996; Bosc et al., 2004).

FIGURE 1

Figure 1. Map of the study area.

Data of occurrence (presence/absence), biomass (metric tons) and abundance (number of individuals) of anchovy and sardine from 2003 to 2016 were retrieved from different sources such as EU-funded MEDIterranean Trawl Survey (MEDITS), EU-funded MEDIterranean International Acoustic Survey (MEDIAS), the Spanish Acoustic Survey “Eco-MEDiterranean” ECOMED and the French Acoustic Survey “Pelagiques MeDiterraneìe” PELMED.

Regarding MEDITS data (Bertrand et al., 2002), species data were collected during demersal trawl surveys carried out during daytime from spring to summer (April to September depending on the country). The months in which the MEDITS survey is carried out (June–July) coincide with the recruitment of sardine, and the spawning activity in anchovy (Palomera et al., 2007). The MEDITS project uses a stratified sampling design based on five depth ranges: 10–50, 51–100, 101–200, 201–500, and 501–700 meters, and GSAs. The surveys in GSA06 were performed by the Spanish Oceanographic Institute (Instituto Español de Oceanografía, IEO), while the French Research Institute for Exploration of the Sea (Institut Français de Recherche pour l’Exploitation de la Mer, IFREMER) conducted the surveys in GSA07.

Pelagic acoustic surveys sampling design in the Mediterranean Sea was homogenized in 2009 under the MEDIAS project (Anonymous, 2008). Prior to 2009, sampling design differed per Mediterranean country. In particular, in the Spanish Mediterranean coast (GSA06) from 2003 to 2008 annual acoustic ECOMED surveys were performed by the IEO during November–December. While acoustics prospections were performed during the day, pelagic trawls were deployed during night time to assess the species composition. The months in which the ECOMED survey was carried out (November–December) coincide with the recruitment of anchovy, and the earliest signs of spawning activity in sardine (Olivar et al., 2001). It worth to be mentioned that fishing during the night or the day does not imply significant differences in surveys catches (Machias et al., 2013). The ECOMED surveys were developed with a systematic design of transects perpendicular to the coastline covering the continental shelf approximately between the 30 and 200 m (Bellido et al., 2008). The PELMED surveys were conducted by IFREMER over the continental shelf of the Gulf of Lions during the summer season (most of the times in July) from 2003 to 2008. All cruises were performed using regularly spaced (12 nautical miles) inshore offshore transects. Acoustic sampling was performed with scientific split-beam echo sounders operating at 38 kHz and calibrated following standard techniques (Foote, 1987). Data were recorded at a constant speed of 8 knots. When the echosounder detected sufficiently long fish traces (≥2 nm) or a change in the echotrace characteristics, a pelagic trawl was deployed to assess the species composition (Van Beveren et al., 2014). The minimum sampling depth varied between 10 and 20 m depending on the area (Ospina-Álvarez et al., 2013). The MEDIAS dataset for both species and the two studied areas (GSA06 and 07) was collected through acoustic samplings from 2009 to 2017 performed during daytime each July by means of scientific split-beam echosounders working at 38 kHz (Saraux et al., 2014). As for the PELMED project, pelagic trawls were performed when the echosounder detected sufficiently long fish traces or a change in the echotrace characteristics. Acoustic data were recorded in both areas following the MEDIAS regular sampling design (parallel transects perpendicular to the coastline) at a constant speed of 8–10 nmi^h–1 (Doray et al., 2010).

All these datasets were available through a European Data Call to the SPELMED project (Coll and Bellido, 2018). It worth mentioning that the only the MEDITS dataset was geo-referenced.

Landing datasets, collected for the fishery harbors located in the studied area, were provided by the General Secretariat of Fisheries of the Spanish Ministry of Agriculture, Food and Environment (MAPAMA) in the GSA06, while for the GSA07 the data was available through a European Data Call to the SPELMED project (Coll and Bellido, 2018).

Environmental Data

For the temporal modeling, three different possible predictors were used to model the habitats of adult anchovy and sardine: Sea Surface Temperature (SST in °C), Sea Surface Salinity (SSS in PSU) and Net Primary Productivity (NPP in mg/m³). They were selected as the most commonly used variables to explain the variability of anchovy and sardine dynamics in the Mediterranean Sea (Fernández-Corredor, 2018).

Monthly averages of oceanographic variables covering the entire study period (2003–2016) were extracted from the Copernicus Marine Environment Monitoring Service¹ for both GSA06 and GSA07.

Bathymetry was obtained from the MARSPEC dataset² with a spatial resolution of 1×1 degree.

To model future scenarios, we extracted different IPCC pathway predictions with a spatial resolution of 1×1 degree from Bio-ORACLE³. In particular, we used the Representative Concentration Pathway (RCP) 2.6 and 8.5 for 2050 and 2100 for SSS and SST.

The RCP2.6 projection scenarios assume least change, with an increase in temperature of 1°C degree in 2050 and 2°C degree in 2100 and of an increase of 0.5 PSU and 1 PSU units for salinity. On the contrary, the RCP8.5 scenarios presume most severe change with an increase in temperature of 1°C degree in 2050 and almost 3°C degrees in 2100, and of 1 and 1.5 PSU units for the salinity, respectively.

Temporal Trends and Correlation Analysis

As preliminary analysis with the attempt to capture the general patterns in the temporal trends of biomass, abundance and landings of the two studied species, we used the “ggplot” package (Wickham, 2011) of the R software (R Core Team, 2019) to plot these indexes. Additionally, in order to assess the complementarity and difference among these different analyzed indexes collected from different oceanographic surveys, we performed a Spearman correlation using the “corrplot” R-package (Wei et al., 2017).

Temporal Modeling

Random Forests (RFs) were applied to assess which environmental variables had influenced the indexes of biomass, abundance and landings of the two studied species historically. We run RFs using data from the different oceanographic surveys (ECOMED, PELMED, MEDIAS, and MEDITS) and landings for both areas (GSA06 and GSA07), using both total values of aforementioned indexes and with the application of a Gaussian distribution. Biomass and abundance data were log-transformed before being introduced in RFs to ensure a normal distribution, as well as to reduce the weight of extreme values. In addition, we tested 12 temporal lags (from January to December) in independent model to assess if anchovy and sardine were significantly influenced by environmental conditions of the months previous to surveys.

RFs are particularly suitable for this task as they do neither require any prior assumptions about the data, nor is their classification accuracy affected by correlations or interactions between variables (Vilela and Bellido, 2015). RFs involve fitting an ensemble of regression trees, first selecting many bootstrap samples from the data, each of which contains ∼63% of the original observations (Strobl et al., 2009). Observations that are not selected in each bootstrap sample are referred to as out-of-bag observations. A regression tree is then fitted to each bootstrap sample, but only a subset of randomly selected variables is used at each node. The trees are fully grown without pruning, and then each tree is used to predict the out-of-bag observations. The out-of-bag estimates are considered a cross-validation of the accuracy of estimates because they are not used in the fitting of trees. The relative importance of each predictor variable is then determined from the misclassification rate for the out-of-bag observations. A sufficiently large number of trees (10,000) were used, and different random seeds were applied to ensure stability in variable importance. The mean squared error (MSE) statistic was computed to measure variables importance (Breiman, 2001). RFs were performed using the “dismo” R-package (Hijmans et al., 2017).

Species Distribution Models

Boosted Regression Trees (BRTs) were used to model the spatial habitat of anchovy and sardine, as this technique deals with non-linear and non-monotonic relationships between response and explanatory variables. BRTs combine regression trees and boosting methods to fit complex non-linear relationships between predictors and a response variable (Elith et al., 2008). BRT parameters were selected to optimize the model using the “caret” (Kuhn, 2008) and “dismo” packages of the R software. The optimal number of boosting trees was assessed with the “gbm.step” function, while the tree complexity of the model was fixed at 2 and the learning rate at 0.01. To account for the spatial autocorrelation in the data, we implemented the residuals autocovariate (RAC) approach (Crase et al., 2012). Spatial autocorrelation was included by adding another term to the model, which represents the influence of neighbor observations on the response variable at a particular location (Escalle et al., 2016; Barcala et al., 2020). For each model, the RAC approach was implemented as follows: first, the model was computed with the selected environmental variables; second, residuals from the selected model were calculated for each grid cell and were used to compute the spatial component by a focal calculation. This allowed cells from a selected neighbor to have a weight of 1 and all other cells a weight of 0. Finally, the spatial component was considered as an explanatory variable in the previous model (Crase et al., 2012). Spatial autocorrelation was tested for each model by calculation of Moran’s index and the Moran statistical test (R package “spdep,” see Bivand et al., 2015), which indicates a correlation between observations depending on the distance between them.

A binomial distribution was implemented to model the response variables of occurrence of the MEDITS dataset. We used these response variables as the only available geo-referenced variables, as well as the large stand most consistent time-series for both areas and species. In addition, temporal trends obtained from MEDITS dataset showed significant correlations with several other available data (ECOMED, PELMED, and MEDIAS; see below).

Validation of the binomial models was conducted through an internal 10-fold cross-validation in which the relationship between occurrence data and environmental variables was modeled using a training dataset (created by a random selection of 75% of the data) and the quality of predictions was then assessed using a validation dataset (created by a random selection of 25% of the data), as advised by Fielding and Bell (1997). The model performance was assessed using the Pearson’s r coefficient, which measures the correlation between predicted and observed values. This coefficient can vary from −1 to 1, with 1 representing a perfect positive correlation between the 2 data sets.

Four future scenarios (RCP2.6 and RCP8.5 for 2050 and 2100) were computed for each species and area (GSA06 and GSA07). Finally, in order to identify potential climate refuges, areas with a probability of occurrence higher than 0.70 (where there is the bulk of occurrence of the species) were selected and mapped.

Results

Temporal Trends and Correlation Analysis

Temporal trends of sardine biomass in the acoustic datasets showed a negative trend for both GSA06 and GSA07 (Figure 2), which was significant in case of GSA07 (rho = −0.65) (Figures 2, 3). Sardine abundance showed positive but non-significant increase in GSA06 and a negative but non-significant decline in GSA07 (Figures 2, 3). Trends of abundance and biomass for sardine were positively correlated in both areas (Figure 3), but significant only in GSA06 (rho = 0.63, Figure 3).

FIGURE 2

Figure 2. Biomass (A–D), abundance (E–H), and landings (I,J) temporal trends from acoustic (ECOMED, PELMED, MEDIAS) and demersal (MEDITS) surveys and official countries landings for European sardine (Sardina pilchardus) in the Geographical Sub-Area 06 and 07 from 2003 to 2016.

FIGURE 3

Figure 3. Spearman’s correlation for European sardine (Sardina pilchardus) among biomass and abundance collected during acoustic (ECOMED, PELMED and MEDIAS) and demersal (MEDITS) surveys, as well as official countries landings, in the Geographical Sub-Area 06 and 07 from 2003 to 2016.

Temporal trends of sardine abundance and biomass using the MEDITS dataset showed a negative decrease in GSA06, with a significant correlation for the biomass (rho = −0.66) (Figures 2, 3). In GSA07 the biomass trend presented a negative trend while the abundance a positive one, but neither trend was significant (Figures 2, 3). Trends of abundance and biomass for sardine in GSA06 and GSA07 were positively significantly correlated (rho = 0.93 and 0.69, respectively; Figure 3).

The declines in landings were large and significant in both areas, with a Spearman correlation value of −0.77 in GSA06 and −0.76 in GSA06 (Figures 2, 3).

Overall, when comparing the three datasets used to describe the temporal trends of sardine in GSA06 and GSA07 (acoustic surveys, trawl surveys and landings), using biomass data from acoustic surveys and from bottom trawling, we observed that in all cases and for both areas sardine showed decreasing trends (Figures 2, 3). These trends were significant for landings in both areas, and for acoustic surveys in GSA07 and trawl surveys in GSA06. It is also interesting to notice that, in all cases, time series were positively correlated. This correlation was significant between landings and acoustic surveys in GSA07 (with rho = 0.50) and between landings and trawl surveys in GSA06 (with rho = 0.53). Between GSAs, we found a significant and positive correlation between landings of sardine (with rho = 0.83) (Figure 3).

Temporal trends of anchovy biomass using the acoustic datasets showed positive increases (Figure 4), significant in the case of GSA06 (rho = 0.84) (Figure 5). Anchovy abundance showed a positive and significant increase with time in GSA06 and GSA07 (rho = 0.71 and rho = 0.61, respectively, Figures 4, 5). Temporal trends of anchovy abundance using the MEDITS dataset showed a positive significant increase in GSA06 (rho = 0.73) and a positive but non-significant decrease in GSA07 (Figure 5). The trends of biomass showed a positive increase in both areas (Figure 4), which was significant in case of GSA06 (rho = 0.78) (Figure 5). Trends of abundance and biomass for anchovy in GSA06 and GSA07 were positively significantly correlated (Figure 5).

FIGURE 4

Figure 4. Biomass (A–D), abundance (E–H), and landings (I,J) temporal trends from acoustic (ECOMED, PELMED, MEDIAS) and demersal (MEDITS) surveys and official countries landings for European anchovy (Engraulis encrasicolus) in the Geographical Sub-Area 06 and 07 from 2003 to 2016.

FIGURE 5

Figure 5. Spearman’s correlation for European anchovy (Engraulis encrasicolus) among biomass and abundance collected during acoustic (ECOMED, PELMED, and MEDIAS) and demersal (MEDITS) surveys, as well as official countries landings, in the Geographical Sub-Area 06 and 07 from 2003 to 2016.

When comparing the three datasets used to describe the temporal trends of anchovy in GSA06 and GSA07 (acoustic surveys, trawl surveys and landings) we observed that, except for landings in GSA07, anchovy showed increasing trends (Figures 4, 5). These trends were significant for landings, acoustics and trawl surveys in GSA06 (with rho’s higher than 0.7). It is also interesting to notice that in several cases time series were correlated. This correlation was highly significant between acoustic surveys and trawl surveys in GSA06 (rho = 0.86), between acoustic surveys and landings in GSA06 (rho = 0.85) and between landings and trawl surveys in GSA06 (rho = 0.81).

Predictors and Lags of Temporal Trends

Summarizing results of temporal modeling with RFs for sardine using biomass, abundance and landings as dependent variables (Figure 6 and Supplementary Figures S1–S3), we observed that NPP showed mainly a negative effect in GSA06 and mixed effects in GSA07. With the exception of the PELMED abundance for which the relationship found with NPP corresponded to the month of the survey, all the relevant relationships found with NPP were with temporal lags of 9–11 months. SST showed a negative impact in GSA07 and mixed effects in GSA06, while SSS showed mainly a negative effect in GSA06 and GSA07.

FIGURE 6

Figure 6. Summary of the Random Forests results for European sardine (Sardina pilchardus) temporal index trends in the Geographical Sub-Area 06 and 07 from 2003 to 2016. This figure shows the ecological and statistical weight of each variable. (+) symbol indicates a positive relationship of the variable with the biomass, abundance or landings and (-) symbol indicates a negative relationship. The colors represent the range, from most to least important variable within the model (red the most relevant, orange medium relevancy and yellow the least relevant) based on the per cent increase in mean squared error- MSE).

Summarizing results of the temporal modeling with RF for anchovy using biomass, abundance and landings as dependent variables (Figure 7 and Supplementary Figures S1–S3), we observed that NPP showed mainly a negative relationship in GSA07, while in GSA06 was less relevant, showing mixed effects. SST showed a positive effect in GSA06 and GSA07, while SSS showed mainly a positive effect in GSA06 and a negative effect in GSA07, especially for landings. In general, results highlighted that time lags were important to explain the relationships between dependent variables for anchovy and sardine, and environmental factors (Figures 6, 7 and Supplementary Figures S1–S3).

FIGURE 7

Figure 7. Summary of the Random Forests results for European anchovy (Engraulis encrasicolus) temporal index trends in the Geographical Sub-Area 06 and 07 from 2003 to 2016. This figure shows the ecological and statistical weight of each variable. (+) symbol indicates a positive relationship of the variable with the biomass, abundance or landings and (-) symbol indicates a negative relationship. The colors represent the range, from most to least important variable within the model (red the most relevant, orange medium relevancy and yellow the least relevant) based on the per cent increase in mean squared error- MSE).

Historical Spatial-Temporal Trends of Anchovy and Sardine Occurrence (2003–2016)

For sardine, RAC-BRTs total deviance explained ranged from 51 to 74% in the GSA07 and from 52 to 78% in GSA06 depending of the year (Supplementary Table S1). Similarly, for anchovy the total deviance explained ranged from 33 to 75% in the GSA07 and from 48 to 78% in GSA06 (Supplementary Table S2). An absence of spatial autocorrelation was detected in all models (Supplementary Tables S1, S2), highlighting that it was accurately handled by the RAC method.

For model validation, reasonably high values for Pearson’s rho were obtained for both species and areas (Supplementary Tables S1, S2). In particular, from RAC-BRTs for sardine in GSA07, rho values between 0.64 and 0.82 were obtained in the cross-validation (Supplementary Table S1). In GSA06 rho values were higher ranging from 0.69 and 0.93. For anchovy in GSA07, rho values ranged from 0.66 to 0.90, while in GSA06 from 0.74 to 0.93 (Supplementary Table S2).

Spatial maps of the occurrence in GSA07 highlighted the preference of sardine for the coastal area and changes with time, with compression of the occurrence niche in 2006, 2012, and 2015 (Supplementary Figure S4). The inter-annual deviation of sardine occurrence from 2003 to 2016 evidenced the large observed changes during the time period, at times more evident in coastal areas or on the shelf and upper slope (Supplementary Figure S5).

In GSA06, spatial maps of the occurrence highlighted a concentration of sardine around the Delta Ebro area, with some years (e.g., 2009, 2013) where the occurrence was especially restricted (Supplementary Figure S6), followed by a range expansion in particular towards the south. The inter-annual deviation of sardine occurrence from 2003 to 2016 also evidenced large changes in the whole study area (Supplementary Figure S7).

For anchovy, RAC-BRTs spatial maps showed that occurrences were mainly located in the coastal and continental shelf area in GSA07 (Supplementary Figure S8). During years 2005, 2006, 2008, 2012, and 2013 the occurrence was closest to the eastern side. Inter-annual variability was also high during the time series (Supplementary Figure S9). In GSA06, spatial maps of the occurrence highlighted an extension of both species’ distribution ranges from 2011 to the present (Supplementary Figure S10). The inter-annual deviation of anchovy occurrence from 2003 to 2016 revealed the expansion of the species in all the study area (Supplementary Figure S11).

Future Spatial-Temporal Scenarios

For sardine in GSA06, projection maps (Figure 8) highlighted a contraction of distribution ranges with only two remaining hot-spots by 2050. These two hot-spots were located (1) north of the Delta Ebro, and (2) south of Cabo de la Nao, in the Gulf of Alicante. These hot-spots emerged for both the optimistic and the pessimistic IPCC pathway scenarios. The hot-spot areas were projected to contract even further in 2100 compared to 2050 predictions. In addition, the map of the variability from the present distribution (Supplementary Figure S12) highlighted that this decrement of the sardine distribution will be gradual from 2050 to 2100 and will occur across GSA06 except for these two hot-spots.

FIGURE 8

Figure 8. Predicted probability of occurrence of the European sardine (Sardina pilchardus) in Geographical Sub-Area 06 for the present situation and using IPCC pathway scenarios obtained with RAC-BRT models.

When areas with a probability of occurrence higher than 0.70 were identified for each scenario, the waters located in the southern area of the Delta Ebro influence resulted as the only zone in which the species may persist in the future (Figure 9). It is worth mentioning that this area will suffer a contraction of the distribution of sardine in the southern part of the Delta Ebro from now to 2100 (Figure 9). This result is very similar when considering the optimistic or the pessimistic IPCC scenario.

FIGURE 9

Figure 9. Climate refuges, hot-spots of the probability of occurrence (>0.7) of the European sardine (Sardina pilchardus) in Geographical Sub-Area 06 for the present situation and using IPCC pathway scenarios obtained with RAC-BRT models. Red color indicates probability of occurrence >0.7, while blue <0.7.

Future projection maps for sardine in GSA07 showed that its spatial distribution and biomass will undergo a gradual decrease from the present to 2100 (Figure 10), with the exception of remaining biomass concentrations in the waters off the Rhone River (Figure 10). The map of the variability from the present distribution to the future scenarios highlighted that this decrease of sardine distribution and biomass will be constant (Supplementary Figure S13).

FIGURE 10

Figure 10. Predicted probability of occurrence of the European sardine (Sardina pilchardus) in Geographical Sub-Area 07 for the present situation and using IPCC pathway scenarios obtained with RAC-BRT models.

When the areas with a probability of occurrence higher than 0.70 where computed in each future scenario, sardine will only persist in the zone around the Rhone River by 2100 (Figure 11). As for sardine in GSA06, the identified area with higher probability of occurrence is very similar in 2050 and 2100, but with a reduction in size in 2100. Results are similar when considering the optimistic or the pessimistic IPCC scenario.

FIGURE 11

Figure 11. Climate refuges, hot-spots of the probability of occurrence (>0.7) of the European sardine (Sardina pilchardus) in Geographical Sub-Area 07 for the present situation and using IPCC pathway scenarios obtained with RAC-BRT models. Red color indicates probability of occurrence >0.7, while blue <0.7.

Anchovy distribution in GSA06 will undergo a decrease in its spatial distribution from the present to 2100. As for sardine, distribution ranges will gradually contract to a minimum area in 2100 (Figure 12 and Supplementary Figure S14). For this species, too, the two main hot-spots that will keep the bulk of biomass will be the area around the Delta Ebro and the Gulf of Alicante. Differently from sardine, these two areas were both identified as hot-spots with a probability of occurrence higher than 0.70 and will persistent into the year 2100 (Figure 13). Results were similar when considering the optimistic and the pessimistic pathway scenarios.

FIGURE 12

Figure 12. Predicted probability of occurrence of the European anchovy (Engraulis encrasicolus) in Geographical Sub-Area 06 for the present situation and using IPCC scenarios obtained with RAC-BRT models.

FIGURE 13

Figure 13. Climate refuges, hot-spots of the probability of occurrence (>0.7) of the European anchovy (Engraulis encrasicolus) in Geographical Sub-Area 06 for the present situation and using IPCC pathway scenarios obtained with RAC-BRT models. Red color indicates probability of occurrence >0.7, while blue <0.7.

Future prediction maps of anchovy in the GSA07 showed very similar results to sardine in the same area probably due to the fact the main changes will be happen all in 2050. For pathway scenarios, prediction and variation maps indicate a restriction in spatial distributions (Figure 14 and Supplementary Figure S15), where the bulk of biomass will only be retained near the mouth of the Rhone River. The area with a probability of occurrence higher than 0.70 will be close to the Rhone river plume and will persist in 2100 although strongly reduced (Figure 15).

FIGURE 14

Figure 14. Predicted probability of occurrence of the European anchovy (Engraulis encrasicolus) in Geographical Sub-Area 07 for the present situation and using IPCC pathway scenarios obtained with RAC-BRT models.

FIGURE 15

Figure 15. Climate refuges, hot-spots of the probability of occurrence (>0.7) of the European anchovy (Engraulis encrasicolus) in Geographical Sub-Area 07 for the present situation and using IPCC pathway scenarios obtained with RAC-BRT models. Red color indicates probability of occurrence >0.7, while blue <0.7.

Discussion

In this study, we present main results of a multi-modeling approach to investigate the main effects of environmental factors on anchovy and sardine biomass, abundance, occurrence and landings. To our knowledge, this is one of the most comprehensive studies about the topic in the Mediterranean Sea.

Regarding changes in sardine abundance and biomass, our results demonstrate a general decreasing temporal trend of biomass of sardine in both GSA06 and GSA07, which coincides with a significant decline in landings in both areas, likely reinforced by the depletion of the commercial stock (Quattrocchi et al., 2016; Coll and Bellido, 2018; Coll et al., 2019). The abundance and biomass of anchovy showed a general increase over time in both GSA06 and GSA07, which coincides with a significant increase of catches in GSA06. On the contrary, we observe a decline in landings in GSA07, which can be motivated by changes in the fishing strategy that shifted from pelagic to demersal target species (Van Beveren et al., 2016), and also due to the fact that most of the pelagic vessels targeting sardine stopped their activity at the end of 2009 due to the stock decline (Fisheries in the Gulf of Lions, UNEP-MAP-RAC/SPA, 2013).

Correlations between trends in acoustic surveys, trawl surveys and landings suggest that these data can be complementary. Overall, correlations between datasets are higher for anchovy than for sardine. These correlations coincide with results from the attempt to compare MEDITS temporal estimates with commercial catches of sardine and anchovy in GSA09 (Ligurian Sea and Northern Tyrrhenian Sea) (Sbrana et al., 2010). For instance, Sbrana et al. (2010) showed that MEDITS survey could be a promising descriptor of anchovy abundance in GSA09 as it is the main targeted species of the fisheries sector in the area. However, it is worth to be mentioned that bottom-trawling is not the best to sampling collector for pelagic species and thus, although similarities were found, results should be taken into account carefully.

Results from the temporal modeling showed mixed effect between variables for anchovy and sardine and primary productivity, especially in GSA06. In the same area, Bellido et al. (2008) found a relatively weak but positive relationship between primary productivity indicators (i.e., Chl-a) and small pelagic species occurrence, suggesting that fish were not directly related to this environmental factor, but that areas with higher primary productivity were more favorable for both sardine and anchovy. However, other studies pointed out that sardine juveniles prefer areas with moderate primary productivity values (Giannoulaki et al., 2011; Tugores et al., 2011) and negative relationships were also found between anchovy adults and primary productivity in other Mediterranean areas (Giannoulaki et al., 2008, 2013; Martín et al., 2008; Quattrocchi et al., 2016). Recent studies point out that primary productivity concentration is negatively correlated with anchovy individuals of age −1 (Bacha et al., 2010; Basilone et al., 2017). Overall, areas with high primary productivity concentration could affect water transparency, increasing the difficulty of finding prey (Fernández-Corredor, 2018). Also, these mixed effects could be due to the different behavior of life stages of the studied species (i.e., juveniles and adults) (Fernández-Corredor, 2018) and could be explained by the high plasticity of their feeding behavior, switching between the selective prey feeding and the non selective filter-feeding mode (Palomera et al., 2007; Zarrad et al., 2008; Ganias, 2009; Costalago et al., 2015). In addition, another possible explication of these mixed temporal effects could due to the fact that the population structure of these spcies has been changing with time. Indeed, in 2006 there were more age classes than in 2016 (Albo-Puigserver et al., 2018).

In GSA06, mixed effects were found among temporal patterns of sardine indexes and temperature, while in GSA07 the relationship was mainly negative. These results are in accordance with literature in which the effect of temperature on sardine indexes was found mainly negative throughout the Mediterranean basin (Bellido et al., 2008; Martín et al., 2008; Katara et al., 2011; Tugores et al., 2011; Martín et al., 2012; Giannoulaki et al., 2013; Brosset et al., 2015; Bonanno et al., 2016; Quattrocchi et al., 2016). Cooler temperature can be indicative of nutrient enrichment processes such as wind mixing, upwelling and river run-off, associated with favorable conditions for fish reproduction and recruitment (Sabatés et al., 2007; Fernández-Corredor, 2018).

On the contrary, for anchovy the relationship with temperature was positive in both areas. This result is consistent with a comprehensive literature review (Fernández-Corredor, 2018). Temperature seems to affect the size of anchovy at the end of the first year of life, and was found positive related with anchovy at age-1 (Bacha et al., 2010; Basilone et al., 2017), and with anchovy body condition factors (Basilone et al., 2006).

Salinity showed overall a negative relationship for both species and areas, with exception for anchovy in GSA06 for which a positive relationship was found with several indexes and temporal lags, as well as with landings. This result is in line with Carpi et al. (2015), which found that marked salinity gradients could have a positive effect on anchovy catches in the Adriatic Sea, while Palomera et al. (2007) show that anchovy is the only small pelagic fish in the NW Mediterranean Sea that spawns in a wide salinity range. On the contrary, landings of sardine in both areas, and for anchovy in GSA07, presented a consistent negative relationship with salinity. Sardine aggregations tend to be higher in shallower areas characterized by lower salinitythat probably are less favorable for this species (Quattrocchi and Maynou, 2017). The negative relationship between salinity and anchovy seems to be correlated to the spawning success (Zorica et al., 2013), lower salinity values endorse the euryhaline condition of anchovies, indicating a connection between salinity fluctuations and their spawning.

Results from the species distribution modeling showed contrasting results between anchovy and sardine and between GSA06 and 07, historically. Anchovy mainly showed an expansion of the population with time in GSA06, especially evident from 2009 onwards, while in GSA07 results illustrate changes between coastal and shelf areas, with some lower than average years in 2005–2008.

Interestingly, in both GSAs for anchovy there was a contraction in the distribution in 2005–2006, coinciding with the years preceding the collapse observed in abundance, landings, growth and body conditions (Van Beveren et al., 2014; Brosset et al., 2015). Also, in both GSAs, in 2013 there was a high retraction of the distribution, coinciding with the year before the second peak of the decline in the species body condition observed 2014 (Albo-Puigserver et al., 2018).

Sardine showed a decrease of occurrence in GSA06 with time, with a slight expansion to the south in the recent years, while in GSA07 a general decline is also observed with a change in concentration of the species from coastal to deeper waters. The last years of historical time series showed a more profound decline. Overall, range contractions during periods of decline in biomass have been observed in a number of different species in the northwest Atlantic including yellowtail flounder (Simpson and Walsh, 2004) and Atlantic cod (Atkinson et al., 1997). The range contraction and expansion are consistent with McCall’s basin hypothesis (MacCall, 1990), where during periods of low biomass, species concentrated in preferred habitats as density-dependent effects declined, and would be in line with recent studies in the Cantabrian Sea (Cabrero et al., 2019). Interestingly, spatially modeled acoustic biomass surveys data for both species in GSA07 obtained similar results (Saraux et al., 2014). This is an indicator that occurrence data from MEDITS could be a complementary descriptor for the distribution of these species.

Also, as for the case of the anchovy, the contraction in the distribution observed in 2013 in both GSAs, coincided with the year before the minimum values of body condition for sardine 2014 (Brosset et al., 2017).

In general, the geographical distribution of the studied species was highly variable from year-to-year, which was expected as these species have a highly mobile behavior and are affected mostly by dynamic environmental factors, which change from year to year. In this sense, averaged environmental functional responses obtained with RAC-BRTs for the entire time series coincide to identify a negative relationship between depth for both anchovy and sardine. These results were in line with several studies that linked anchovy and sardines’ indexes to shallower waters (Giannoulaki et al., 2008; Carpi et al., 2015; Bonanno et al., 2016), showing preferences for areas shallower than 100 m (Bellido et al., 2008). Tugores et al. (2011) found a higher probability of occurrence for sardine in waters less than 65 m depth during summer in GSA06. In GSA07, Saraux et al. (2014) identified preferential areas close to the coast for anchovies (isobaths of ~70 to 100 m) and near the coast and in the Western part of the Gulf of Lions for sardines. RAC-BRTs also captured a mainly positive relationship with primary productivity, while the relationships with temperature and salinity were not linear and showed some thresholds in agreement to what has been found in the literature (Palomera et al., 2007; Fernandez et al., submitted).

Interestingly, and despite the different historical trajectories of anchovy and sardine in the study area, future pathway projection for these two species highlighted that both anchovy and sardine will undergo a reduction in their spatial distributions in both GSA06 and GSA07 due to future climate conditions. Due to the mid trophic level position of sardines and anchovies and their key roles in the ecosystem of the western Mediterranean Sea (Coll et al., 2008, 2019), a reduction on their populations could impact the ecosystem dynamics across entire marine food web, from species that mainly prey on these forage fish to the dynamics of plankton (Tudela and Palomera, 1999; Navarro et al., 2017; Coll et al., 2019). For both species, optimal environmental values (SST and SSS) will be exceeded, first for sardine but eventually for anchovy, too, in combination with a reduction of eggs and larvae of anchovy in the last decades in GSA06 (Maynou et al., 2020) since the reproduction success of sardine and anchovy are highly influenced by SST and SSS condition (Palomera et al., 2007). The areas that will keep a probability of occurrence higher than 0.70 are around the Rhone River in the GSA07 and around the Ebro river in the GSA06 for both sardine and anchovy, and waters around the Gulf of Alicante for anchovy. This coastal area close to the river delta, might be preferred areas in the future, due to the river runoff that could compensate the changes expected from an increase of SST, with increasing stratification of the water column affecting the plankton dynamics to a lower biomass or changes in composition (Calvo et al., 2011). Optimistic and pessimistic IPCC pathway scenarios result in very similar future outcomes for both species and this is probably due to the fact that the most drastic changes will already occur in 2050 for these species and will continue until 2100.

Nowadays pathway scenarios did not include how the river waters itself will change in the future. Future studies should assess whether and how river waters and land use will affect species distribution predictions. However, the identified areas could be seen as “future climate refuges” for these species. Identifying areas that are most likely to persist under adverse climate change conditions in the future is important for anticipating the complex impacts that climate change will have on ecosystems and societies (Monsarrat et al., 2019). The socio-ecological implications of some species that recolonize part of their historical niche or occupy new areas should be carefully considered, in particular for species with a high socio-ecological interest such as small pelagic fishes.

These results are especially relevant to elaborate spatially explicit management plan (e.g., Marine Spatial Planning, MPA, etc.) for small pelagic species. Indeed, they show that despite the high mobility of small pelagic fish and high interannual variability, some areas consistently offer favorable habitats and should be prioritize for conservation measures and marine spatial planning.

Data Availability Statement

The data analyzed in this study is subject to the following licenses/restrictions: The dataset was received from an European Data Call for the development of SPELMED project and we are not authorized to share it. Requests to access these datasets should be directed to “European Data Call” or MC, bWNvbGxAaWNtLmNzaWMuZXM=, and JB, am9zZW0uYmVsbGlkb0BpZW8uZXM=.

Author Contributions

MP and MC designed the research. MP, EF-C, and MC performed the research. MP, EF-C, and JS analyzed the data. All authors wrote the manuscript.

Funding

This study was carried out within the Spanish Research project PELWEB (CTM2017-88939-R) funded by the Spanish Ministry of Science, Innovation and Universities, the European Research Contract SPELMED (EASME/EMFF/2016/032) funded by the ECEASME and the Catalonian Government PELCAT projects (CAT 152CAT00013 and TAIS ARP059/19/00005). The information and views set out in this publication are those of the author(s) and do not necessarily reflect the official opinion of the Executive Agency for Small- and Medium-sized Enterprises (EASME) or of the European Commission. Neither EASME, nor the Commission can guarantee the accuracy of the data included in this publication. Neither EASME, nor the Commission or any person acting on the EASME’s or on the Commission’s behalf may be held responsible for the use which may be made of the information contained therein.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

Thanks are due to the General Secretariat of Fisheries of the Spanish Ministry of Agriculture, Food and Environment (MAPAMA) for providing fishery landings data from GSA06. Thanks are also due to European institutions to provide GSA07 fishery landings data and survey data throughout a European Data Call for the development of SPELMED project.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2020.00622/full#supplementary-material

Footnotes

References

Albo-Puigserver, M., Navarro, J., Coll, M., Layman, C. A., and Palomera, I. (2016). Trophic structure of pelagic species in the northwestern Mediterranean Sea. J. Sea Res. 117, 27–35. doi: 10.1016/j.seares.2016.09.003