Fisheries-independent data come from MEDITS research surveys, carried out annually during May-September along the northern Mediterranean, using an experimental bottom trawl to assess demersal fishing resources (Bertrand et al., 2002; Spedicato et al., 2019). The sampling strategy of these research surveys follows standard protocols (MEDITS group, 2017). Data used in the present study covered the period 2001–2020, the depth range 40-800 m and the Geographic Sub-Areas (GSAs) 1, 2, 5 and 6 ( Figure 2A ), established by the General Fisheries Commission for the Mediterranean (GFCM), which encompass the westernmost range of the Mediterranean basin (Iberian Peninsula and Balearic Islands). More information about MEDITS sampling methodology is included in Supplementary Information.

Data on standardized density (individuals/km²) of deep-water rose shrimp were compiled from a total of 3,822 sampling hauls (1,025 species presences and 2,797 absences; Supplementary Table S1 ). These samples were aggregated in a 0.1x0.1° grid (ca. 100 km²) and the proportion of presence and abundance of the species per each cell and year was calculated, to avoid bias derived from varying sampling effort over time. Subsequently linear regression with time (year) was computed per each cell and only significant positive slopes (p-value<0.05) were retrieved. The calculations were done using MATLAB v.R2010b and plotted with “M Map: A mapping package for MATLAB”, v.1.4m (www.eoas.ubc.ca/~rich/map.html).

As fisheries data, landings of deep-water rose shrimp along the period 2001–2021 were used to estimate nominal and standardized Catch Per Unit of Effort (CPUE), following the methodology by Soto et al. (2002) as a proxy of its abundance. More detailed information about the methodology applied to estimate CPUE is included in Supplementary Information . The nominal CPUE has variations, mainly associated with changes in the characteristics and composition of the fishing fleet (factors related to catchability; Hilborn and Walters, 1992) and environmental factors (factors related to abundance; Sund et al., 1981). To avoid biases related to these aforementioned factors, catch indices were standardized (Beverton and Holt, 1959) by means of Generalized Linear Models (GLM; Punt et al., 2000). The GLM included the nominal CPUE as a response variable, being the explanatory variables factors associated with abundance, catchability and random error.

Fishing landings were analysed for the most important fishing harbours along the whole Mediterranean coast of Iberian Peninsula, grouped by sub-regions ( Figure 2A ): Catalonia (Roses and Tarragona), Gulf of Valencia (GoV: Cullera and Denia), Gulf of Alicante (GoA: Santa Pola and Villajoyosa) and Andalucía (Almería and Velez-Málaga). Data from the Fish Auction of Palma (Balearic Islands), where all the fishing landings from Mallorca harbours are commercialized, were also used. Nominal CPUE was fitted using a gamma distribution with a logistic link and explained deviations of 46.3%. Finally, the nominal CPUE was replaced by the CPUE estimated by the model, obtaining the standardized CPUE.

Nine environmental variables were considered as potential predictors ( Table 1 ), including main previously suggested drivers of habitat preference for deep-water rose shrimp: temperature and salinity (Sbrana et al., 2006; Fortibuoni et al., 2010) and surface chlorophyll, as proxy of biomass of primary producers in the water column (Guijarro et al., 2009). More information about these variables is included in Supplementary Information.

Data were mostly extracted from Copernicus Marine Environment Monitoring Services (CMEMS; https://marine.copernicus.eu/), which use remote sensing observations and include reanalysis products specifics for the Mediterranean ( Table 1 ). Monthly data were averaged for the spring-summer season, matching the timing of MEDITS surveys (May-September) and all variables were reprocessed to match the same grid of 0.1x0.1° spatial resolution used for biological data (see Section 2.1). All environmental variables were explored for collinearity, outliers and missing data, before being used in the analysis and modelling (Zuur et al., 2009). Finally, values of the selected environmental variables were extracted from all MEDITS hauls locations where the species was captured, using extract function from “raster” R package (Hijmans, 2023).

A set of simulations (ensemble) were gathered from four downscaled Regional Climate Models (RCMs) ( Supplementary Table S2 ), to describe potential future conditions in the western Mediterranean, under the frame of the Coupled Model Intercomparison Project Phase 5 (CMIP5). This sub-basin presents complex oceanographic features (e.g. fronts, deep water formation) that need to be properly captured in the models, through using high-resolution simulations (Ruti et al., 2016; Schoeman et al., 2023).

We followed previous similar efforts (e.g. CERES project; Maynou et al., 2020) for the Representative Concentration Pathways (RCPs) 4.5 and 8.5 greenhouse gas concentration scenarios. RCP4.5 and RCP8.5 represent moderate and high greenhouse gas concentration, that increase radiative forcing on Earth by 4.5 and 8.5 Wm⁻², respectively, at the end of the 21st century. All the simulations outputs were interpolated to the same grid (0.1x0.1°) and bathymetry. The climatological values corresponding to May-September were extracted and averaged for 2006–2020 (present, representing the current conditions), middle 21st century (2040–2050, coded as 2050s) and the end of the 21^st century (2090–2100, coded as 2100s). We followed similar widely implemented approaches for projecting shifting spatial distribution of marine species (Assis et al., 2018).

Bias corrected maps were computed subtracting delta values per cell from original projection maps. Delta values are the difference between the original projection of 2006–2020 (“present”) and the same overlapping period from the physical reanalysis of the Mediterranean ( Table 1 and Supplementary Table S2 ). This reanalysis included data assimilation and properly captured the main hydrological and biogeochemical features of this area (Ramirez-Romero et al., 2020). Furthermore, the performance of each ensemble member ( Supplementary Table S2 ) was assessed using the physical reanalysis as reference and collected in a Taylor diagram (Taylor, 2001) ( Supplementary Figure S4 ). All the model details and validation could be found in supplementary information.

To avoid overfitting and to reduce the correlation and collinearity between variables, we filtered the initial set of predictors based on two metrics. First, all predictor covariates were examined using Pearson’s rank correlation (Wood, 2006). Pairs of variables with high correlation values (Pearson correlation; r>0.7) were identified and only one of the correlated pairs was included in the modelling process (Lezama-Ochoa et al., 2017). Subsequently, multicollinearity between predictor variables was tested by calculating the variance inflation factor (VIF) with a cutoff value of 5 (Hahlbeck et al., 2017; Lopez et al., 2020) and using the corvif function from the “AED” package in R (Zuur et al., 2009). As a result, highly correlated and multicollinearity variables were tested separately during the modelling, and only the one with the highest total explained variance was included in the final models.

For the spatial distribution modelling we used a two-stage approach to predict habitat suitability in density of deep-water rose shrimp, as the survey data was inflated to zero (Barry and Welsh, 2002). Generalized Additive Models (GAM; Hastie and Tibshirani, 1990; Wood, 2006) were used to build two models: (i) a binomial model (presence/absence) to predict the occurrence of the species, which was fitted through a binomial distribution with a logistic link function; and (ii) an abundance model (truncated; all zeros excluded), which was fitted through a Gamma distribution with a logarithmic link function.

The GAM models were fitted in the R software using the “mgcv” package (Wood and Wood, 2023), where the degrees of freedom of the smoothing functions were limited to 10 for the latitude and longitude interaction spline (spatial component) in the binomial model and for the abundance model, achieving the best fits. For the explanatory variables, the degrees of freedom of the smoothing functions were limited to 5 in bathymetry and 4 in sea bottom temperature. A forward step-wise variable selection procedure and the best model was selected based on Akaike Information Criterion (AIC), the unbiased risk estimator (UBRE) (Burnham and Anderson, 2004), significant explanatory variables and the highest explained deviance. These two models were combined by multiplying the predictions from both steps to obtain the final predicted value or preferred habitat model.

The evaluation of the predictive model of density was carried out in two phases: (i) for the model of occurrence and abundance, correlating the occurrence and the abundance predicted and observed by using the Pearson correlation coefficient (r); and (ii) for the binomial model by using the area under the receiver operating curve (AUC) (Fielding and Bell, 1997; Elith et al., 2006), and specificity, sensitivity, and True Skill Statistic (TSS) (Allouche et al., 2006). This model validation was performed using the “PresenceAbsence” R package (Freeman and Moisen, 2008).

The spatial structure of the forecasts, derived from the occurrence models in the three time windows (present, 2050s and 2100s) and taking into account the two different warming scenarios (RCP4.5 and RCP8.5), were compared using Schoener and similarity statistics (Warren et al., 2008). These analyses were performed using the “NicheOverlap” function of the “dismo” package (Hijmans et al., 2011) within the R software. In addition, we compared abundance and density maps using Pearson’s r spatial correlations with the “corLocal” function (“R raster” package; Hijmans, 2023). The spatial distribution model, evaluation current predictions and validation projections details could be found in Supplementary Information.

The analysis of MEDITS data for the period 2001–2021 showed that the initial presence of the deep-water rose shrimp was mainly restricted to the southernmost part of the study area (GSAs 1, 2 and southern part of GSA 6; Figure 3A ). However, density increasing trends were registered along the two-decades period at the central and northernmost areas (GSA 6) and noticeably around the Balearic Islands (GSA 5; Figure 3B ). The CPUEs of this species shared the same patterns. Initially, the highest catches were found in the southernmost basin and since 2007 a rising trend was registered in central and northernmost areas (GSAs 5 and 6), reaching similar levels in the whole area during the recent years ( Figure 2B ).

According to the cross-correlation matrix of the environmental variables included in the models and their variance inflation factors: Sea Bottom Temperature (SBT) was highly correlated with Sea Bottom Salinity (SBS), Integrated Phytoplankton Carbon Biomass and Integrated Chlorophyll were also highly correlated, bathymetry was highly correlated with SBS ( Supplementary Figure S1A ) and VIF values indicated multicollinearity between them ( Supplementary Figure S1B ).

Bathymetry, SBT and spatial component (latitude-longitude interaction) were the main drivers of the occurrence and abundance of the deep-water rose shrimp, with the greatest explained deviations (46.3 and 50%, respectively; Table 2 ). The most important covariate to explain the occurrence and abundance of the species was latitude-longitude interaction, with an explained deviation of 32.7 and 47.9%, respectively, followed by bathymetry with 29.8 and 19.3%, respectively and finally SBT with 12.8 and 17.2%, respectively.

For the predictive performance of species occurrence, the Pearson correlation coefficient (r) ranged from 0.58 to 0.59 and for abundance it was 0.40-0.41 ( Supplementary Table S1 ). In addition, the occurrence model indicates a good degree of discrimination between presences and absences (AUC>0.5), good ability to fit true positive and negative predictions (SPEC and SENS ranged between 0.5 and 0.8) and good degree of similarity between the occurrence of the species and the available evidence (TSS ranging between 0.2 and 0.55). The spatial autocorrelation of the residuals of the final models indicates the absence of significant spatial autocorrelation in the residuals of the models of presence (Moran’s I= 0.11; p-value= 0.45) and abundance (Moran’s I= 0.69; p-value= 0.24). The species had occurrence and abundance maxima between 200 and 500 m depth ( Figure 4 ). Regarding the SBT, both GAM models showed positive linear responses, increasing probabilities with higher SBT, particularly from ca. 15.5°C. Predictive maps showed an increasing northward gradient along the north-western Mediterranean, mainly due to SBT. All predictive maps for the density of the species are shown in Supplementary Figure S2.

SBT patterns matched the deep-water rose shrimp distribution shift ( Figure 3 ; Supplementary Figure S3 ). During the last decades (2001–2020), the southernmost area (GSA 1) was the warmest, but the central and northern areas (GSAs 5 and 6) presented the highest linear warming trends ( Supplementary Figure S2 ).

Overall, the ensemble members (four RCMs) reproduced the observed patterns at present conditions (2006–2020) reasonably well, particularly for spatial correlations ( Supplementary Figure S4 ). The ensemble averaged maps (bias-corrected) showed progressive warming of SBT along all the study area in both scenarios ( Supplementary Figures S5 , S6 ), more marked in the most extreme scenario (RCP8.5; Supplementary Figure S6 ). The internal variability in the ensemble (captured by the standard deviation) was generally larger in the 2050s compared to the 2100s ( Supplementary Figures S5 , S6 ).

Deep-water rose shrimp are expected to slightly increase its occupied area, with an expansion from 50 to 60% of the occupied area. Furthermore, projected densities showed a noticeably increase for both RCPs, rising from current estimated densities (10⁴ individuals Km^-2) to 10⁵ individuals Km^-2 (2050s) and reaching 10⁶ individuals Km^-2 at the end of the 21st century (2100s, Supplementary Figures S5 , S6 ). This density raise was projected in almost the entire study area, but it was more marked in the central and northernmost parts (GSAs 5 and 6). Figure 5

Simulations registered a good level of agreement regarding the internal coherence of the ensemble, excepting a more divergent response under RCP4.5 in 2050s. Although being the most temporally remote period, in 2100s both RCP presented low internal variability within the ensemble ( Figure 6 ).

Regarding the robustness of these forecasts, in general the similarity statistics highlighted a good level of overlap in the occurrence predictions of the two data sets. In fact, Schoener’s D and Warren’s I ranged from approximately 0.87 to 0.99 for the RCP4.5 and RCP8.5 scenarios ( Supplementary Table S3 ). For the abundance and density maps, the spatial correlations also indicated a high similarity between the predictions of different models in of present, 2050s and 2100s for the RCP4.5 and RCP8.5 scenarios ( Supplementary Figures S7 , S8 ).

Ensemble of simulations presented a high degree of match regarding the spatial forecast of the occupied area and density of the species ( Supplementary Figure S9 ). In addition, the increase in density is also spatially forecasted with a reasonably internal coherence, particularly for the furthest future conditions and for the RCP8.5 scenario ( Supplementary Figure S10 ).

The Mediterranean Sea, a climate change hot spot (Shaltout and Omstedt, 2014; Tuel and Eltahir, 2020), presents unique oceanographic processes typical from marginal seas, such as the thermohaline circulation and deep-convective water formation, tightly coupled to intense atmospheric processes (Millot and Taupier-Letage, 2005). Because of this, a high sensitivity to climate forcing is expected in this basin, as small changes in the water masses hydrology can induce large variations in the thermohaline circulation and energy or mass budgets. In addition, complex mesoscale dynamics (e.g. fronts and filaments) are broadly present over the basin (Mason et al., 2023).

Despite this complexity, RCMs could properly capture these main features and coordinated multi-model and multi-scenario studies were carried out in an initiative as Med-CORDEX, specifically for the Mediterranean (https://www.medcordex.eu/; Ruti et al., 2016). According to previous RCM ensembles, which provide information on the range of projection outcomes, the water temperature in the intermediate layers (150-600 m depth) of the western Mediterranean, where many of its key demersal fishing resources inhabit, are expected to increase in 1.32 ± 0.16 and 2.23 ± 0.42°C for RCP4.5 and RCP8.5, respectively (Soto-Navarro et al., 2020). Recent works even forecasted a collapse of the deep-water formation in the Gulf of Lions by 2040–2050, attributed to an increasing stratification potentially leading on pervasive modifications of the basin circulation (Parras-Berrocal et al., 2022).

The positive linear response of the deep-water rose shrimp to SBT ( Figures 4B, D ), modeled in this study from research survey data, confirmed the fisheries-data presented here ( Figure 2 ) and the results by Mingote et al. (2024) in the northernmost area using landings. This species has been favored in terms of abundance and density, due to the water warming affecting the Mediterranean during recent decades (Colloca et al., 2014) and, therefore, can be considered a sentinel of climate change (Hazen et al., 2019). Given that the species is an important component of demersal crustacean communities (e.g. Abelló et al., 2002) and fishing resources (Bertrand et al., 2002; Spedicato et al., 2019), this finding is particularly useful for early detecting warming hotspots in deep ecosystems, where tracking changes in communities is challenging.

Although water temperature above seabed has been validated as a good descriptor of habitat suitability for this species in the western Mediterranean (Abelló et al., 2002; Colloca et al., 2014; Mingote et al., 2024), our forecasting results could be considered in the future for their improvement. The spatial distribution of demersal species can be projected more accurately and precisely by incorporating mechanistic/physiological or density-dependent effects, among other factors. Moreover, the predictions were carried out without taking into account other factors such as fishing effort, biotic interactions (e.g. competition and predation, reproduction/fecundity rates) and ecological processes, such as dispersal limitations or habitat connectivity. This simplification may limit the accuracy of predictions, especially for ecologically complex systems.

The warming progression along the Mediterranean Sea can extensively alter the availability of fishing resources (e.g. Ben Lamine et al., 2022), emerging as an essential factor for the assessment and management of fisheries. It becomes especially relevant in areas like the Mediterranean, where most demersal stocks have been overexploited for long time (Vasilakopoulos et al., 2014), generally impaired by a high fishing pressure, combined with the low selectivity of bottom trawling (Colloca et al., 2013). As a consequence, most commercial species are currently overexploited (FAO Fishery and Aquaculture Global Statistics, 2023) and show truncated size- and age-structures of their populations, which can cause then the loss of resilience in the face of climate change (Hilborn and Walters, 1992; Pauly and Zeller, 2016; Barange et al., 2018). In fact, the synergistic effects of fishing-induced demographic changes and climate variation on Mediterranean fish population dynamics have been already described (Hidalgo et al., 2011).

Understanding medium and long-term changes in marine ecosystems allows fishing stakeholders to anticipate potential changes in exploited resources in the coming decades. This enables them to propose adaptation strategies to minimize climate change impacts and regulate fishing efforts. These strategies are crucial, as climate-driven changes in the distribution, abundance and/or migration of species might not align with fishing efforts, due to the complex response of the socio-economic system (Barange et al., 2018). This adaptation, defined as a “process of adjustment in ecological, social or economic systems to actual or expected climate and its effects” (IPCC, 2014), has been insufficiently considered for Mediterranean fisheries. In fact, National Adaptation Plans regarding climate change generally show limited references to adaptation measures for fisheries (Hidalgo et al., 2018).

The predictions and projections on the impact of climate on fisheries must be made at regional and local scales, where historical data series can contribute to improving complex models. These scales are particularly suitable for adapting fisheries industries and fisheries management, since global models are of limited value (Brander et al., 2013). By focusing on regional and local scales, it is possible to tailor strategies to the specific conditions and challenges of different areas. This allows for more precise and effective adaptation measures that take into account local ecosystem dynamics, socio-economic factors and community needs. Adapting at these scales ensures that fishing sector can respond more accurately to the impacts of climate change at different regions, leading to more sustainable and resilient fisheries. In this context, the present results are even more crucial considering that the deep-water rose shrimp is one of the six target species under Regulation (EU) 2019/1022 of the European Parliament and the Council, which aims to sustainably manage demersal stocks in the western Mediterranean. The implementation of this multi-annual plan is leading to a reduction up to 40% in the fishing days of the bottom trawl fleet during the recent years. According to our results, this reduction would not be fully justified in the case of this species.

In conclusion, this study improves the understanding of the habitat preferences and current and future distribution of one of the most relevant demersal species in the western Mediterranean. Our results indicate a northward expansion and an increase in the density of the species, mainly driven by warming trends of water above seafloor. These findings are fundamental to develop effective fisheries management strategies in the face of global warming, ensuring the sustainability of this key fishing resource. The use of a RCM ensemble has proven to be advantageous in providing robust and ecologically interpretable projections, ultimately contributing to better scientific advice for fisheries planning adapted to environmental changes caused by global warming.