Resilience dynamics and productivity‐driven shifts in the marine communities of the Western Mediterranean Sea

Abstract Ecological resilience has become a conceptual cornerstone bridging ecological processes to conservation needs. Global change is increasingly associated with local changes in environmental conditions that can cause abrupt ecosystem reorganizations attending to system‐specific resilience fluctuations with time (i.e. resilience dynamics). Here we assess resilience dynamics associated with climate‐driven ecosystems transitions, expressed as changes in the relevant contribution of species with different life‐history strategies, in two benthopelagic systems. We analysed data from 1994 to 2019 coming from a scientific bottom trawl survey in two environmentally contrasting ecosystems in the Western Mediterranean Sea—Northern Spain and Alboran Sea. Benthopelagic species were categorized according to their life‐history strategies (opportunistic, periodic and equilibrium), ecosystem functions and habitats. We implemented an Integrated Resilience Assessment (IRA) to elucidate the response mechanism of the studied ecosystems to several candidate environmental stressors and quantify the ecosystems’ resilience. We demonstrate that both ecosystems responded discontinuously to changes in chlorophyll‐a concentration more than any other stressor. The response in Northern Spain indicated a more overarching regime shift than in the Alboran Sea. Opportunistic fish were unfavoured in both ecosystems in the recent periods, while invertebrate species of short life cycle were generally favoured, particularly benthic species in the Alboran Sea. The study illustrates that the resilience dynamics of the two ecosystems were mostly associated with fluctuating productivity, but subtle and long‐term effects from sea warming and fishing reduction were also discernible. Such dynamics are typical of systems with wide environmental gradient such as the Northern Spain, as well as systems with highly hydrodynamic and of biogeographical complexity such as the Alboran Sea. We stress that management should become more adaptive by utilizing the knowledge on the systems’ productivity thresholds and underlying shifts to help anticipate both short‐term/less predictable events and long‐term/expected effects of climate change.


| INTRODUC TI ON
Resilience has become a conceptual cornerstone in marine management and a popular narrative for marine conservation as it provides the opportunity to communicate ecosystem health and the capacity to overcome both natural and anthropogenic disturbances (Darling & Côté, 2018). However, its applications and implications are still partially disconnected from the ecological basis (Johnson & Lidström, 2018). This is largely due to the difficulty to find metrics and methods to apply efficiently the resilience theory in an empirical context (e.g. Collie et al., 2004;Scheffer, 2009but see, Vasilakopoulos et al., 2017Sguotti et al., 2019;Capdevila et al., 2020).
Resilience is generally defined in two complementary ways in ecological literature: ecological and engineering resilience. Ecological resilience is the capacity of a system to absorb disturbance and reorganize while undergoing change so as to retain essentially the same function, structure, identity and feedbacks (Walker et al., 2004). By contrast, engineering resilience is an expression of the time required for an ecosystem to return to an equilibrium or steady state following a perturbation (Holling, 1996). While engineering resilience focuses on the stability around a single equilibrium, ecological resilience emphasizes conditions far from the equilibrium, where a system may flip to an alternate domain (Holling, 1996). In any case, both the ability of ecosystems to resist and absorb disturbance without changing basin of attraction (ecological resilience) and their rate of recovery towards a prevailing steady state (engineering resilience) are complementary emergent properties to be secured for the oceans' sustainability.
Herein, 'resilience' refers to ecological resilience, unless otherwise specified, as it is based on the theory of critical transitions and alternate stable states (Scheffer, 2009). Complex natural systems with eroded resilience, ranging from populations to ecosystems and socioecological systems, can exhibit abrupt reorganizations in their components in response to external stressors (e.g. Litzow & Hunsicker, 2016;Scheffer, 2009;Scheffer et al., 2001). These system reorganizations are broadly known as regime shifts, defined as large and abrupt changes in the structure and function of complex systems. Different types of shifts can be identified attending to the relationship between the ecosystem state(s) and the fluctuations of the driver(s), broadly categorized as continuous or discontinuous (Scheffer et al., 2001;Selkoe et al., 2015).
Regime shifts alter the properties of ecosystems by triggering persistent changes in their structure, dynamics and internal feedbacks, which can affect specific ecosystem functions (Möllmann et al., 2015). To operationalize resilience in ecosystem-based management, it is crucial to identify which functions, and associated species, are (un)favoured in alternate regimes. Recent studies have suggested, for instance, that pelagic and thermophilic species of smaller size and of low trophic level will benefit from sea warming scenarios in temperate ecosystems (e.g. Moullec et al., 2019). While large-scale thermal gradients shape geographical variation in marine traits (Pecuchet et al., 2017), regional variation is adapted to local environmental variability. Identifying which traits and functions are favoured in the alternative stable states before and after a regime shift would allow to identify specific management measures tailored to the new ecosystem and environmental realities at regional scale.
This may help to avoid chronic and profound impacts that could hamper ecosystem services, particularly in environmentally heterogeneous systems.
Beyond ocean warming, climate change impacts are also associated with less predictable changes in environmental conditions, as well as with extreme weather events such as heatwaves or strong winter storms. While extreme events are mainly known to cause mass mortality episodes, they can also cause a reorganization of the ecosystem components (Maxwell et al., 2019). However, in ecosystems with low resilience, small changes in the environmental conditions can also trigger drastic changes in the prevailing regime with substantial modification in key functions and the most frequent traits in the communities. Such environmental drivers include also changes in primary production (i.e. productivity shifts) that can cascade fast in the food web and trigger strong modifications on the upper trophic levels (Fu et al., 2018).
The Mediterranean Sea is also one of the most exposed areas worldwide to the numerous impacts of climate change, including one of the highest warming rates, and increasing frequency of extreme weather events, and also subjected to productivity changes (Macias et al., 2018). It is also exposed to multiple anthropogenic impacts, such as overfishing and invasive species (e.g. Coll et al., 2012). The cumulative effects of different impacts and the patched structure of marine populations and communities make Mediterranean systems highly responsive to environmental fluctuations, and their resilient capacities highly driven by climate variability including primary production changes (Piroddi et al., 2017).
Here, we hypothesize that putative recent climate-driven shifts are occurring in the Western Mediterranean associated with changes in the relative contribution of different groups of species more adaptive by utilizing the knowledge on the systems' productivity thresholds and underlying shifts to help anticipate both short-term/less predictable events and long-term/expected effects of climate change.  Vasilakopoulos & Marshall, 2015;Vasilakopoulos et al., 2017) to assess the response type of the ecosystems against several candidate environmental stressors, quantify ecological resilience and its temporal variation (i.e. resilience dynamics), and build stability landscapes with alternate valleys of attraction.

| Systems and data
We focus on two environmentally contrasting demersal ecosystems: in Northern Spain and in the Northern Alboran Sea (hereafter referred as 'Alboran'; Figure 1). The Northern Spain system corresponds to the eastern Balearic Sea associated with the Iberian Peninsula coast and is characterized by large-scale environmental gradients, including surface and intermediate water masses temperature, and primary productivity regimes resulting from the confluence of deep ocean convection from the Gulf of Lions, riverine nutrient inputs by the Ebro river and sub-mesoscale and frontal processes (Ramirez-Romero et al., 2020). The Alboran is, by contrast, a transitional system between the Atlantic and the Mediterranean ecosystems with strong hydrodynamism characterized by two permanent anticyclonic gyres encircled by the strong eastward flow and a higher mean primary productivity compared to Northern Spain.
This study used information obtained from a scientific bottom trawl survey (MEDITS program; note that ethical approval was not required for this research), conducted to evaluate the status of demersal resources in the early summer (May-June) of every year from 1994 to 2019. Here, we focused on the demersal and megaepibenthic community of the continental shelf (<200 m; Figure 1), having sampled 1,436 and 391 hauls in the Northern Spain and the Alboran, respectively, with an average depth of 89 m and 87 m. Data used in this study were retrieved at species level, averaging biomass per swept area (i.e. kg/km 2 ) caught for each species and sampling station in every year.
With regards to environmental data, we focused on three potential drivers of the biological system: chlorophyll concentration (Chl a), sea surface temperature (SST) and an integrated Regional Hydroclimatic Index (RHI). These variables were selected attending to observed species-specific responses and low trophic level effects in previous studies in the same areas, in both pelagic and demersal ecosystems (e.g. Puerta et al., 2014). We retrieved seasonal records of environmental variables to calculate two seasonal averages: for spring (March-May) and winter (December-February), assuming that species more sensitive to changes in primary production may be more associated with the winter signal as it is associated with the main phytoplankton bloom. SST information at 1/16 degree resolution for 1992-2019 was obtained from the Mediterranean physics re-analyses available at the Copernicus Marine Service (https:// marine.coper nicus.eu).
To calculate a Regional Hydroclimatic Index (RHI), several environmental variables were obtained from the NCEP-DOE Reanalysis 2 fields provided by NOAA (https://psl.noaa.gov/data/gridd ed/data. ncep.reana lysis2.html) for 1992-2019. Anomaly fields were retrieved at a monthly scale for air temperature, sea surface temperature, atmospheric sea level pressure, 500 hPa geopotential height and precipitation rates. The RHI quantifies an integrated hydroclimatic variability at regional scale synthesizing the aforementioned variables by means of the first axis of a Principal Component Analyses (PCA). High RHI values are associated with higher atmospheric sea level pressure and 500 hPa geopotential height, and negative values with high precipitation rate, both in Northern Spain and in Alboran.

| Integrated Resilience Assessment
The IRA is a three-step methodological framework which applies the concepts of resilience and folded stability landscapes in an empirical multivariate context through the combination of multivariate analysis, non-additive modelling and a resilience assessment (Vasilakopoulos et al., 2017). This way, the IRA elucidates the system dynamics and shift mechanisms in response to external stressors.
As a first step, a PCA was applied to the time-species matrix of each studied area to identify the main modes of community variability. The PCAs were based on the correlation matrices of the species' biomasses, following log-transformation. PCA converts a set of possibly correlated variables into new linearly uncorrelated composite variables (principal components-PCs) using an orthogonal transformation of the data to a new coordinate system (Vasilakopoulos & Marshall, 2015). The standardized value that each data point has in the PC coordinate system (PC-score) is then calculated. As the majority of individual species reported in the two ecosystems in 1994-2019 were caught only sporadically, we applied the PCA on six datasets for each area, including species caught in all 26 years or at least in 25, 24, 23, 22 or 21 years, respectively. This way, it was examined how the perception of the system dynamics would change in response to the number of species analysed while ensuring that only well-sampled species would be retained. That 'optimal' PCA selected was the most parsimonious one (i.e. with the lowest number of species), as long as it captured similar dynamics with the ones including more species.
The sequential regime shift detection method (STARS), modified to account for temporal autocorrelation (Rodionov, 2006), was applied to detect significant shifts in the mean values of PC1 and PC2 of the optimal PCAs. STARS estimates a Regime Shift Index (RSI), that is, a cumulative sum of normalized anomalies relative to each value of the time series analysed, and uses it to test the hypothesis of a regime shift occurring in that year (Rodionov, 2006). Here, we used a cut-off length of 3 years and a significant probability thresh- The composite variable (PC1 or PC2) retained to be further analysed within the IRA was selected based on its capacity to capture directional changes and non-overlapping states, and on the degree of its coupling with the environmental variables. To identify the key species driving the community trends through time, the contributions of all species on PC-scores (loadings) of the retained PCs were checked.
To investigate the type of the relationship between the biological system, as captured by PC1 or PC2, and each of the environmental stressors (winter and spring Chl a, SST and RHI), the fits of relevant generalized additive models (GAMs) and non-additive threshold GAMs (TGAMs) at 0-to 2-time lags were compared. The effect of the stressors on the system was examined at 0-to 2-year lags to account for a potential delay in the environmental effect on the community sampled given that environmental effects typically influence spawning and early life stages, while the sampled biomass is usually dominated by specimens older than 0 years old. This maximum of 2 years is in line with previous Mediterranean studies (Damalas et al., 2021;Tsimara et al., 2021;Vasilakopoulos et al., 2017), and it was selected because of the narrow demography of Mediterranean fish species (Colloca et al., 2013) and their generally low generation time. For fish, Merluccius merluccius was the species with the highest generation time and above 2 years, with the rest of the species likely below this value (expect elasmobranchs). Crustaceans and cephalopods have life cycles of 2-3 years and 1 year, respectively.
Generalized additive models assume additive and stationary relationships between the response and explanatory variables, while TGAMs can represent an abrupt change in the relationships between the response and explanatory variables (i.e. a threshold) in a specific year supporting the occurrence of a regime shift (Ciannelli et al., 2004). In other words, a GAM describes a system that changes in a continuous way in response to the corresponding change of its stressor(s), while a TGAM represents a system response curve that is folded backwards, forming a fold bifurcation with two tipping points ( Figure S1). Accordingly, we fit 18 GAMs and 18 TGAMs (3 stressors × 2 seasons × 3 lags) for each system, using either PC1 or PC2 as a response variable. To compare the goodness of fit of GAMs and TGAMs and select the optimal model, we computed the 'genuine' cross-validation squared prediction error (genuine CV; gCV), which accounts for the estimation of the threshold line and the estimation of the degrees of freedom for the functions appearing in all additive and non-additive formulations (Ciannelli et al., 2004). Additionally, a sensitivity analysis was carried out using also PC3 as response variable for the GAM-TGAM comparison, to confirm that PC1 or PC2 captured best the system response to environmental stressors in both systems.
In both Northern Spain and Alboran, the branches of the respective optimal TGAMs identified earlier (i.e. the ones with the lowest gCV values) were considered as the attractors forming the fold bifurcations, in accordance with the theory of critical transitions (Scheffer et al., 2001). Subsequently, the resilience of each annual system state within the fold-bifurcation was quantified based on the distance of each state from its attractor and tipping point (Vasilakopoulos & Marshall, 2015). Specifically, to estimate annual resilience values (Res y ), the sum of the horizontal distance of each annual system state from the tipping point of its regime minus the vertical distance of each annual system state from its fitted attractor (TGAM branch) was used ( Figure S1). This calculation requires the x-and the y-axes of the bifurcation diagram to be measured at the same scale; hence, both the environmental stressor and the system PC were standardized (mean of 0 and standard deviation of 1; Vasilakopoulos et al., 2017). To calculate the position of the tipping point of each regime along the trajectory of its respective attractor, the x-coordinates of the tipping points were set so as to ensure that the lowest Res y estimate within each regime was equal to zero. Finally, Res y was scaled by dividing with the maximum value observed to calculate relative resilience (rRes y ). The stability landscape with its alternate basins of attraction emerged through linear interpolation of all rRes y values onto a 100 × 100 grid.

| Species indicators relevant to taxonomic, functional and temperature affinity
To quantify which species groups were favoured ('winners') and un- give an equal weight to all species and avoid a higher influence of species with higher mean biomass, log-transformed biomass of each species was standardized (mean of 0 and standard deviation of 1).
To evaluate the potential increase in species with higher affinity to warmer temperatures, we used the information of the AQUAMAPS database (https://www.aquam aps.org) to extract the optimal temperature of the habitat distribution for each species, using only the area of distribution of species with probability of occurrence higher than 50% and weighing the temperature by the probability of occurrence. For the few species for which the information was not available in AQUAMAPS, we calculated their preferred temperature based on the mean value of the presence data recorded in the Ocean Biodiversity Information Systems (OBIS, https://obis.org). Finally, we calculated the community-weighted mean temperature (CWMT; Cheung et al., 2013) for each year as Regarding community-weighted mean temperature (CWMT) of Northern Spain, a significant increase was observed in the whole time series for species of commercial interest (sl = 0.006, p < 0.05, Figure 6). This trend was consistently observed across all taxonomic groups, while only significant in fish and crustaceans (sl = 0.004 and sl = 0.025, respectively, both p < 0.05, Figure 6).
In Alboran, unlike in Northern Spain, significant long-term decreasing and increasing trends were observed for non-commercial and commercial species, respectively (sl = −0.02 and sl = 0.04, both p < 0.05), mainly referring to sparids and carangids (Figure 4e; Table S2). This is consistent with the community shift observed in the early 2000s that, in respect to fish, was partially related to a decrease in the contribution of opportunistic fish species (Figure 4g).
Some less pronounced change was also observed in equilibrium and periodic fish species. The temporal development of crustaceans and cephalopods reflected the early 2000s shift with a general decrease in pelagic invertebrates (mainly cephalopods) and an increase in both benthic cephalopods and crustaceans ( Figure 5), with a significant F I G U R E 2 Temporal variation of the main modes of the community variability revealed by the yearly PC-scores for the first two principal components of the PCAs applied for the Northern Spain (above) and Alboran (below). Time series of PC-scores of the two modes (a, c) and their biplot (b, d) trend in the latter (sl = 0.004, both p < 0.01, Figure 5). The estimates of the CWMT also showed a general increase for commercial species over the study period (sl = 0.005, p < 0.05, Figure 6), while a sharp decrease associated with the early 2000s shift was observed in the CWMT of the non-commercial species. The increase in CWMT was mainly associated with fish (sl = 0.005, p < 0.05, Figure 6). Noncommercial CWMT decrease was associated with crustacean species ( Figure 6f).

| Integrated Resilience Assessment
PC1 and PC2 captured better than other PCs the directional underlying trends in Northern Spain and Alboran, respectively, and exhibited non-overlapping regimes (Figure 2). Additionally, the examination of all continuous (GAMs) and discontinuous (TGAMs) fits using PC1 or PC2 as a response variable and all candidate stressors at 0-to 2-year lags as explanatory variables, revealed that the a TGAM of PC1 with winter Chl a at a 2-year lag and a TGAM of PC2 with winter Chl a at an 1-year lag provided the best fits for Northern Spain and Alboran, respectively (lowest gCV; Table 1; Figure S6). The sensitivity analysis using PC3 as a response variable in both areas did not yield any stronger coupling with environmental stressors (Table 1). Hence, PC1 was retained as system indicator (PC1sys) for F I G U R E 6 Community-weighted mean temperature (CWMT, °C) of the demersal community in Northern Spain (left) and Alboran (right): commercial (black line) and non-commercial species in the whole community (a, b), fish (c, d), crustaceans (e, f) and cephalopods (g, h). Dotted black lines indicate that the temporal trend is significant over the whole period (p < 0.05; only statistically significant temporal trends are presented). Grey background represents the new ecosystem states after the regime shifts identified, that is, after 2009 and 2001 in Northern Spain and Alboran, respectively

TA B L E 1
The optimal (lowest gCV) GAM or TGAM models identified using PC1, PC2 or PC3 as a response variable and Chl a, SST and RHI of spring and winter at 0-to 2-year lags as explanatory variables in the Northern Spain and Alboran systems. In bold, the model exhibiting the best fit for each area. All PCs were standardized prior to the GAM-TGAM fitting, thus making gCVs comparable F I G U R E 7 Stability landscapes and estimated relative resilience for the ecosystems of Northern Spain (a) and Alboran (b). The folded stability landscape of each ecosystem exhibited two basins of attraction and two tipping points (F1, F2). Continuous black lines indicate the attractors, dotted black lines indicate the possible extensions of the attractors, grey dashed lines indicate the boundaries of the basins of attraction, black and red arrows indicate the pathways of the regime shifts. Note that in (a) no resilience estimation was carried out for 1994, 1995, as they appear to be the 'tail' of an older state, and for 2008 due to its corresponding attractor being unclear In both Norther Spain and Alboran, the current system states appear to be more resilient at lower concentrations of Chl a compared to the states pre-shift.

| DISCUSS ION
The quantification of the dynamics of ecological resilience in two complex and contrasting communities of the Western Mediterranean Sea elucidated unknown ecosystem mechanisms associated with recent shifts. Intriguingly, the ecosystems investigated, while very complex and exposed to a variety of stressors, were both found to be largely driven by fluctuations in oceanic productivity with a temporal lag of 1-2 years. While different patterns were generally observed across taxonomic and functional groups in the two ecosystems, certain consistencies also emerged. Invertebrate species of short life cycle were generally favoured in recent chlorophyll-a concentration regimes, particularly benthic species in the Alboran.
By contrast, opportunistic fish species were relatively unfavoured in the two ecosystems in the recent periods, with an increase in commercial species of short life cycle in the Alboran and equilibrium and periodic species in the Northern Spain. Beyond the productivity dynamics driving demersal community in the two systems, our results also suggest the subtle and long-term effects of sea warming and relative decrease in fishing. Notably, the analysis of survey data allowed a more integrated representation of both resilience and ecosystem dynamics compared to previous IRAs applied to commercial catches (Damalas et al., 2021;Ma et al., 2021;Tsimara et al., 2021;Vasilakopoulos et al., 2017).
Our results suggest that the coupling between the phenology and favourable environmental conditions (e.g. primary production blooms) is critical at the community level. While the two areas studied here are associated with different productivity regimes with a higher spatial complexity of the phytoplankton blooms in Alboran Beyond consistent short-term responses to primary productivity, other long-term drivers likely affect systems synergistically, increasing its sensitivity to environmental variability. Sea warming is the most evident long-term stressor, causing an increase in the biomass of species with higher temperature affinities and/ or reduction of those with lower affinity (e.g. Vasilakopoulos et al., 2017). In our studied ecosystems, the community-weighted mean temperature indicated a similar increasing pattern for the fish community, but this trend was only mirrored by invertebrates (crustaceans and cephalopods) in Northern Spain. Crustaceans and fish species of low commercial interest and high thermal affinities decreased after early 2000s in Alboran, along with the increase in some species of low thermal affinities. This is likely due to high turnover and replacement of species fostered after the early 2000s in the Alboran Sea that often receives the species runoff from the Atlantic side of the Iberian Peninsula (Real et al., 2021). Additionally, the impact of fishing has been suggested as a key chronic stressor in the Mediterranean increasing the ecosystems' sensitivity to productivity changes (Damalas et al., 2021;Fu et al., 2018). However, a recent decrease in fishing effort and displacement of fishing activity to deeper strata in the Western Mediterranean allowed a relative recovery of highly vulnerable species in the shelf such as elasmobranchs (Ramírez-Amaro et al., 2020) and benefitted regional fish diversity (Farriols et al., 2019). This is consistent with our results in Northern Spain and the increase in the relative contribution of some periodic commercial and equilibrium species in recent years, and also with the recent increase in commercial species of short life cycle in the Alboran Sea.
The concurrent long-term warming, and moderate and heterogeneous decrease in fishing impact could be favouring species of short life span in recent chlorophyll-a regimes. This seems to be the case in Northern Spain, where the long-term increase in the relative contribution of species with higher thermal affinity seems to be widespread, as expected in the Mediterranean Sea and particularly for generalist groups (Doubleday et al., 2016;Moullec et al., 2019). With the change in the demersal communities' regimes, species of short life span, such as some benthic cephalopods and crustaceans, were favoured in Northern Spain and particularly in Alboran suggesting a fuelling of the benthopelagic coupling with an increased efficiency of energy transfer to the benthic ecosystems. By contrast, opportunistic fish species were not equally favoured. This could be due to the high fish diversity in the Mediterranean, particularly for species with opportunistic life-history strategy, which represent a high proportion in the Mediterranean fish communities (Pecuchet et al., 2017). These aspects could be of more important in systems of highly hydrodynamic and of biogeographical complexity such as the Alboran Sea (Real et al., 2021). Opportunistic species are likely to suffer from high replacement and turnover rates, as they tend to have specialized behaviour and trophic interactions and could be easily replaced under the combination of warming and changes in the productivity regime. The Alboran Sea also has an additional element of complexity due to its high hydrodynamics and low environmental predictability to which resident species and community assemblages have adapted over their evolutionary history (Winemiller & Rose, 1992), with opportunistic species becoming more sensitive and replaceable under fluctuations of primary productivity. This is consistent with the fact that the first axis of the PCA in the Alboran showed no directional trend, suggesting an efficient niche replacement and reduced impact of species coexistence.
Emerging ecosystem shifts such as the productivity-driven shifts reported here or shifts driven by sporadic strong perturbations history traits composition has also been detected, switching traits systems to shorter life span, smaller maximum size, higher optimal temperature and shallower optimal depth (Tsimara et al., 2021). This new traits configuration could be making communities more sensitive to regional changes in primary productivity.
The regime shifts described in this study comply with the properties defined by the theory of critical transitions in complex ecological systems (Scheffer et al., 2001;Selkoe et al., 2015) and resemble other climate and/or fishing-driven transitions displaying strong hysteresis (Litzow & Hunsicker, 2016;Möllmann et al., 2015).
Nevertheless, seldom such shifts have been associated with productivity fluctuations as in our case. The presence of hysteresis implies that a new ecosystem state remains resilient and fails to reverse to its previous state even after a perturbation is reversed (Beisner et al., 2003;Scheffer et al., 2001). This kind of system dynamics implies a difference in the set of variables controlling the ecosystem response in the alternate states, with a potential change of the main environmental forcing on primary production ( Figure S8). Positive feedbacks such as changes in the trophic relationships could contribute to amplify the amount of change in the system.
Most marine stewardships worldwide are currently working towards climate resiliency in ecosystems management (e.g. Holsman et al., 2019). To do that, a general systematic approach for improving resource management by learning from management outcomes (i.e. adaptive management) has become the means towards resilient socioecological landscapes (Walker et al., 2004). Management approaches which do not properly integrate across spatiotemporal scales are maladapted to unidirectional change and extreme events (Holsman et al., 2019). Our results are in this sense paradigmatic as they reveal region-specific resilience dynamics to short-term unexpected events and also the long-term effects of sea warming and fishing decrease. However, further research is needed to operationalize the results obtained in this study. First, early warning indicators specific for productivity-driven regime shifts need to be developed (Litzow & Hunsicker, 2016). Second, these indicators should be combined with the identification of the specific areas in which fine-scale processes affecting productivity can cascade over a broad spatial and/or temporal extent (Selkoe et al., 2015). Moreover, the productivity thresholds for regime shifts identified here could be compared with forecasts of productivity in the studied areas to facilitate resilient socioecological landscapes under adaptive management frameworks. In other words, future regime shifts could be anticipated by managers to help them adapt their policies accordingly, combining strategies to manage both long-term warming impacts and shortterm events such as productivity fluctuations and shocks.

ACK N OWLED G EM ENTS
We thank members of the ICES working group COMEDA for their input, stimulating discussions and feedback during the preparation of this paper. We also thank L. Pecuchet for providing the database with life-history strategy scores for fish. MEDITS data collection was co-funded by the EU through the European Maritime and Fisheries

DATA AVA I L A B I L I T Y S TAT E M E N T
PCA axes and chlorophyll-a data that support IRA analyses are available at https://doi.org/10.5281/zenodo.5750293 (Hidalgo, 2021).