An indicator-based approach for assessing marine ecosystem resilience

Marine ecosystems are under threat from a range of human pressures, notably climate change, overexploitation, and habitat destruction. The resulting loss of species and biodiversity can cause abrupt and potentially irreversible changes in their structure and functioning. Consequently, maximizing resilience has emerged as a key concept in conservation and management. However, despite a well-developed theory, there is an urgent need for a framework that can quantify key components promoting resilience by accounting for the role of biodiversity. In this study, we applied an indicator-based approach to assess the potential resilience of marine ecosystems using the North Sea as an illustrative case study. More specifically, we quantified and compared multiple indicators of ecological resilience, estimated based on high-resolution monitoring data on marine demersal fish species, combined with information on ecological traits. Our results show a pronounced spatial structuring of indicators, including both similarities and differences among individual metrics and indicators. This implies that high resilience cannot be achieved by maximizing all individual aspects of resilience, simply because there seems to be inherent trade-offs between these components. Our framework is generic and is therefore applicable to other systems and can inform spatial planning and management.


Introduction
The diversity of life on Earth, whether represented by variation in genes, species, or habitats, is currently threatened by multiple human pressures, including climate change, overexploitation, and habitat destruction (Cardinale et al., 2012;Hansen et al., 2016;IPBES, 2019). Consequently, the abundance and spatial distribution of species worldwide are rapidly changing (Barnosky et al., 2011;Pimm et al., 2014;Ceballos et al., 2015;Carmona et al., 2021). The loss of species and biodiversity can cause abrupt and potentially irreversible changes in the structure and functioning of ecosystems and can ultimately lead to alternative ecosystem states (Scheffer et al., 2001;Folke et al., 2004;Bianchelli et al., 2016). Such large-scale changes have been observed over a broad range of ecosystems, including in terrestrial, freshwater, and marine systems (Scheffer et al., 2012;Conversi et al., 2015;Knowlton N. 1992). Based on the empirical support for the occurrence of alternative stable states in nature, ecological resilience has emerged as a highly relevant concept to understand and manage biodiversity loss to ensure the integrity and stability of ecosystems in the face of change (Walker, 1995;Folke et al., 2004;Conversi et al., 2015;Chambers et al., 2019).
Resilience was formally introduced by Holling (1973) as "a measure of the persistence of systems and of their ability to absorb change and disturbance and still maintain the same relationships between populations or state variables". Since then, the term has been increasingly used within the scientific community and society at large, but conflicting and contradictory meanings of the concept have emerged. Consequently, Holling later distinguished between two types of resilience: engineering and ecological resilience (Holling, 1996). Engineering resilience is defined as the return time to equilibrium, or a measure of the rate at which a system approaches equilibrium. It relies on the core assumption that the system in question has a single steady state around which it may fluctuate, but always strives to return to after a perturbation (Folke et al., 2004). Conversely, ecological resilience relies on the assumption that the system in question may contain multiple stable states and is defined as "the amount of disturbance that a system can withstand without changing self-organized processes and structures" (Peterson et al., 1998;Folke et al., 2004;Walker et al., 2004;Chambers et al., 2019). Approaches for quantifying ecological resilience cover a broad range of fields from spatial assessments (Cumming 2011;Gladstone-Gallagher et al., 2019;Gunderson, 2000), functional diversity assessments (Angeler and Allen, 2016), and identification of thresholds and early warning signals (Standish et al., 2014, Burthe et al., 2016. Despite the well-developed theory and concepts of resilience, translating ecological resilience into quantifiable and directly measurable metrics is challenging (Dakos et al., 2012;Scheffer et al., 2012;Angeler and Allen, 2016). This challenge inevitably hampers progress towards fulfilling the ambitious policy goals to protect biodiversity and maintain resilience, functions, and services of ecosystems (Balmford et al., 2010;Steffen et al., 2015;IPBES, 2019;Hermoso et al., 2022). Consequently, there is an urgent need for a framework that can identify and quantify key components promoting ecological resilience by accounting for the role of biodiversity. Figure 1. A conceptual illustration of ecosystem states with high and low resilience and the four components contributing to resilience, as well as the underlying data used to estimate these.
The exploration of biodiversity, notably the diversity of species traits and their links to ecosystem functions, allows for identification of key components of ecological resilience (Norling et al., 2007;Violle et al., 2007;Mouchet et al., 2010;Letters, 2017;Delgado-Baquerizo et al., 2020). It has been demonstrated that biodiversity is related to ecosystem resilience, stability, and functioning (Buchmann, 2006;Cardinale, 2006;Cardinale et al., 2012), especially over large spatio-temporal scales (Loreau et al., 2001;Duffy et al., 2007;Stachowicz et al., 2007;Hillebrand and Matthiessen, 2009). The cross-scale resilience model provides a conceptual framework for understanding how ecological resilience emerges in complex systems (Peterson et al., 1998). The concept does not measure resilience directly but highlights several key biodiversity-related components that potentially promote resilience over time, including the richness, redundancy, evenness, and response diversity of species within and between functional groups (Figure 1). These components represent complementary biodiversity metrics that have been shown to be significant determinants of ecological resilience (Giller et al., 2004;Gagic et al., 2015;Angeler and Allen et al., 2016;Sundstrom et al., 2018). This provides an important framework to explore and evaluate a set of potential indicators of resilience.
Here, we applied an indicator-based approach to assess the potential resilience of ecosystems. Specifically, we provided a proof of concept by quantifying and comparing multiple indicators of ecological resilience in a large marine ecosystem. The set of metrics selected for this study can all be estimated based on available monitoring data on species abundances and traits. In this study, we use the North Sea as an illustrative case study. The North Sea region is heavily impacted by a range of human activities (Halpern et al., 2008;Bowler et al., 2020) and has previously demonstrated substantial declines in numerous marine fish stocks (Fernandes and Cook, 2013;Clausen et al., 2018;Lindegren et al., 2018), as well as abrupt and large-scale changes in the overall structure and functioning of the ecosystem (Beaugrand, 2004;Alheit et al., 2005;Lindegren et al., 2012). This highlights the importance of including aspects of resilience in regional marine and fisheries management and conservation, preferably within a holistic ecosystem-based assessment of the ecosystem (Worm et al., 2009;Dickey-Collas et al., 2010;Rickels et al., 2019).

Survey data
We compiled high-resolution survey data of fish species occurrences and abundances from the North Sea International Bottom-Trawl Survey, made available by the International Council for the Exploration of the Sea (https://datras.ices.dk). We extracted all available samples (i.e. unique hauls) from the first quarter of the year (January-March) from 1983 to 2020 ( Figure S9). For each haul, species were identified, the number of individuals was determined, and their length was measured. To standardize for sampling effort and remove strictly coastal species not adequately sampled by the gear, we only retained hauls with a duration between 25 and 35 min and a depth below 20 m. Furthermore, we restricted the analysis to demersal and other bottom-dwelling fish species, excluding pelagic species. All species names were checked and validated according to the World Register of Marine Species (https: //www.marinespecies.org). As spatial resolution, we used the standard ICES rectangles (1 • longitude by 0.5 • latitude grid cells) commonly used in fisheries management. Hauls from all years were assigned to the corresponding grid cell. The final dataset contained 14470 unique hauls and 161 species distributed across 38 years and 185 rectangles. As there was heterogeneity in the number of hauls between grid cells and years, we created species accumulation curves (SACs) to standardize for differences in sampling effort. For each grid cell, an SAC was calculated using the "specaccum" function in the vegan R package (Oksanen et al., 2013). SACs give an insight into how much of the total species richness in an area is covered by a given level of sampling effort (Gotelli, 2001;Chao et al., 2009). Then, we fit nonlinear Michaelis-Menten curves to each of the SACs and estimated the minimum number of hauls required to reach 50, 65, or 80% of the derived asymptotic species richness for each grid cell. Based on the SAC, we ran 100 randomizations of hauls with a sample size given by the fitted SACs for each grid cell (Keating and Quinn, 1998;Maureaud et al., 2019).

Trait data
In addition to the species abundance data, we used a comprehensive dataset containing trait information for marine fish species recorded in the bottom-trawl survey . The dataset contains information on 14 traits, among which 9 are continuous and 5 are categorical traits. Traits are defined as measurable features affecting the fitness of organisms through the processes of feeding, reproduction, and survival (Violle et al., 2007;Beauchard et al., 2017;Butt and Gallagher, 2018). Traits can help determine a species' response to environmental changes and provide insight into the functional role of the species in the ecosystem (Foden, 2013). Traits can be morphological (size and body shape), physiological (metabolic or growth related), or behavioural (migration and feeding patterns). Each trait and the trait extraction procedure are described in detail in . For the analyses, we selected traits with most available data to keep as many species as possible, which includes the following traits: offspring size, age of maturity, fecundity, length infinity, growth coefficient, maximum length, and maximum age ( Table 1). The trait information was used to assign each species into functional groups sharing a similar set of traits. We adopted a similar approach to Mouillot et al. (2007), where reef fish species were grouped into "functional entities" constituting unique combinations of traits at different levels of resolution. Since, in our case, the traits consist of continuous variables, the functional grouping was achieved through hierarchical clustering of the provided traits for all species included in the analysis based on Gower distances using the FD function in the "dbFD" package (Gower, 2008;Villéger, 2008). To test the sensitivity of results to the level of functional grouping, we per-formed the analysis on three hierarchical levels of functional groups amounting to 10, 20, and 50% of the total number of species. At each grouping level, we measured and compared the set of resilience indicators.

Indicators of resilience
To investigate and compare spatial patterns in the potential resilience of North Sea fish communities, we estimated a range of biodiversity-related indicators and metrics that are considered key indicators of resilience, namely functional group richness (Richness), functional group redundancy (Redundancy), functional group evenness (Evenness), and functional response diversity (Response diversity) ( Table 2). Richness was estimated as the mean number of functional groups present per grid cell, based on the previous classification of species into functional groups according to their traits. Redundancy was estimated as the mean number of species present per functional group. Evenness was represented by two complementary metrics, including the mean of Pielou's evenness index (EvennessPE) as well as the mean of Shannon diversity index (EvennessSh) (Heip, 1974;Spellerberg and Fedor, 2003) per functional group. Response diversity was also represented by two complementary metrics: (i) the mean degree of synchrony/asynchrony in species abundance fluctuations within functional groups (Response diversitysyn) (Gonzalez and Loreau, 2009;Lindegren et al., 2016) and (ii) the Portfolio effect, representing the mean variance within each functional group (Response diversitypoe) (Anderson et al., 2013). To estimate the metrics, we used the R packages "synchrony" (Loreau and Mazancourt, 2008) and the function "synchrony", as well as the package "ecofolio" and the "pe_mv" function (Anderson et al., 2013). Although the potential resilience of a system may change over time, some of the metrics used in this study, notably for response diversity, require long-term time-series to assess the degree of synchrony or asynchrony in species fluctuations. Therefore, our analysis is limited to a spatial comparison among indicators and metrics. Hence, all indicators and their associated metrics were estimated for each grid cell and represented as the mean and standard deviation across all years based on 100 random bootstraps of hauls with a sample size given by the fitted SACs for each grid cell.

Community resilience
To compare and investigate how the different indicators of resilience co-vary throughout the study area, we performed a principal component analysis (PCA) using the resilience indicators for each grid cell as input. Subsequently, we investigated and visualized the degree of correlation among indicators in the first two dimensions. Additionally, we mapped the associated loadings (correlations) of each grid cell in space for PC1 and PC2. Furthermore, indicators calculated at each of the three levels of functional groupings were included in the PCA to assess the sensitivity of results to the level of functional grouping.
To represent a potential metric of overall community resilience, we calculated a "joint resilience indicator" as the sum of standardized z-scores (i.e. zero mean and unit variance) across all individual indicators and mapped it out in space. The calculation follows a standard approach for constructing composite indicators commonly used in sustainability assessments (OECD 2008; Singh et al., 2009; Tarasewiczs et

Species trait Unit Description
Offspring size mm Offspring size represents the egg diameter for fish, length of egg case for skates and rays or body length of a new-born pup for sharks Age at maturity year Age at which 50% of the population becomes mature Fecundity No. eggs /year Number of eggs or offspring produced per year by female Length infinity cm Based on the von Bertalanffy model (Pauly et al., 1987), which creates the growth curve for a time series used to model mean species length from its age Growth coefficient 1/year von Bertalanffy growth coefficient K Maximum length cm Maximum species length recorded Maximum age years Maximum species age recorded Where:  Higher values for each metric denote potential higher resilience and vice versa, with the exception of Syn, where higher values (i.e. close to 1) indicate lower response diversity and therefore lower potential resilience. Please note that the output (shown in the result section) represents mean values across functional groups per grid cell after 100 random bootstraps. al., 2021) and serves to highlight areas of generally higher or lower overall scores. Although the approach allows for different weighting factors for each individual indicator, we calculated an unweighted index. This is because we have no prior knowledge about the degree to which different indicators and metrics used in our study reflect the true, yet unmeasured, resilience of the system. However, to give equal weighting and balance the number of variables, we chose to include one metric for each indicator when computing the joint resilience indicator. This was achieved by identifying and removing the highly correlated, and therefore redundant, PE and Poe metrics for evenness and response diversity, respectively, based on the preceding PCA. Finally, to identify potential regional differences in resilience indicators, we calculated a distance matrix based on the mean value of each resilience indicator. Then, we used the function "chclust" from the R package "rioja" to identify potential clusters within the study area. To highlight differences in the individual resilience indicators between clusters, we created boxplots for each metric.

Spatial patterns of resilience indicators
The geographical patterns of the resilience indicators all show large spatial variability throughout the study area (Figure 2, Supplementary Figure S9). Richness shows a higher number of functional groups (>12) present primarily in the northern part of the North Sea, Skagerrak, and Kattegat. Relatively few functional groups (<10) are present in the central part ( Figure  2a). In terms of redundancy, areas with low redundancy (<2) are found in the central northern part of the North Sea. Higher values are found in the southern part, Kattegat, and around the United Kingdom. This indicates a higher number of species per functional group in the latter areas. The highest redundancy values were found around the English Channel (>3) (Figure 2b). In general, there is a low evenness throughout the area (Figure 2c and d), which indicates that the abundances between species within functional groups are unevenly distributed. The lowest evenness scores (EvennessSh < 0.2 and EvennessPE < 0.1) are located at the centre of the North Sea around 55-58 • N and in the Skagerrak and Kattegat. Evenness increases towards the South and around the UK coast (EvennessSh > 0.3 and EvennessPE > 0.12). The two metrics of response diversity (response diversitysyn and response diversitypoe) show considerably more heterogeneity throughout the area compared to the other indicators (Figure 2e and f). In terms of response diversitysyn, lower values (<0.6) are generally found in the mid/northern part of the North Sea, indicating that this area has a higher degree of asynchrony in species fluctuation within functional groups. Conversely, response diveristysyn increases towards southern areas, and the highest values (>0.8) are found in the Kattegat (Figure 2e). Likewise, response diversitypoe is heterogeneously distributed and does not show any clear overall spatial patterns, except a tendency towards slightly higher number of grid cells with high values in the central part of the area (Figure 2f). The overall spatial patterns of the resilience indicators are similar and insensitive to the three hierarchical levels of functional groupings (Supplementary Figures S1 and S2). However, we note that even if the spatial patterns are similar, then the range of values differs for some metrics, notably for richness where the level of functional grouping determines the overall number of functional groups present.

Community resilience patterns
The first two principal components of the PCA explain 30.64% (PC1) and 23.08% (PC2) of the total variation among all indicators. Evenness and redundancy load similarly on PC1 as they have high and positive loadings on PC1 ( Figure  3). Conversely, richness and response diversity are grouped together and demonstrate high loadings on PC2. This demonstrates a clear distinction between indicators of functional redundancy and evenness and those representing functional group richness and response diversity. As evident from the PCA biplot, there is a high degree of correlation between complementary metrics for each resilience indicator, as well as between metrics calculated based on each of the three levels of functional groupings. As an example, both metrics of evenness (EvennessPE and EvennessSh) load similarly on PC1, regardless of functional grouping. Furthermore, there is a negative correlation between response diversitypoe and response diversitysyn, reflecting that the metrics are inversely related, i.e. representing asynchrony and synchrony, respectively. To investigate the spatial covariation of the indicators of resilience, we plotted the loadings of PC1 and PC2 in space. For PC1, the highest loadings are found towards the southern part of the North Sea and the western part of the area along the UK coast (Figure 3b). High PC1 values indicate that most of the variation within the communities can be explained by evenness and redundancy, which load strongly on PC1, whereas low-value areas are driven by response diversity. For PC2, the values are higher in the northern parts of the study area and Kattegat, while the lowest values are found in the central North Sea (Figure 3c). The high values on PC2 indicate that most of the variation within the communities is due to richness and response diversitysyn, which load positively on PC1, whereas low-value areas are associated with response diversitypoe. Similar overall spatial patterns in resilience indicators are also shown using the joint resilience indicator that demonstrates high values (>5) primarily in the South/West North Sea, Skagerrak, and Kattegat ( Figure 4). The lowest values are found in the central and northern parts. There are no marked changes in the overall spatial patterns using the different levels of functional grouping. Finally, the cluster analysis based on the resilience indicators shows four significant spatial clusters ( Figure 5). The first cluster (red) is located in the southern part of the North Sea, while the second cluster (blue) is situated primarily in the Skagerrak/Kattegat area. The third cluster (green) is in the northwestern area, while the fourth (purple) is in the northeastern part of the area. In terms of their underlying indicators, the first area can be characterized by high redundancy, evenness, and response diversitysyn ( Figure 6). The second and third clusters are somewhat similar, but with a slightly lower evenness and a slightly higher redundancy and response diversitysyn in the latter area. The fourth area has the lowest richness, redundancy, evenness, and response diversitysyn overall. The distribution of resilience indicators is broadly consistent across all levels of functional groupings except for the 10% cutoff, which demonstrated a different separation of clusters in the southern North Sea (Figure 5a).

Discussion
As biodiversity loss accelerates worldwide (Pimm et al., 1995;Butchart et al., 2010;Cardinale et al., 2012), it is important to improve our understanding of how biodiversity enhances ecological resilience and ensures the integrity and stability of ecosystems in the face of global change. Despite a solid theoretical understanding of ecological resilience (Folke et al., 2004;Holling, 1973) and readily defined metrics to measure its underlying components (Peterson et al., 1998;Angeler and Allen, 2016), no formal empirical quantification and comparison among the various metrics has been undertaken to characterize the potential resilience in large marine ecosystems (Scheffer et al., 2015;Ingrisch and Bahn, 2018;Chambers et al., 2019). This inevitably hampers our progress towards fulfilling ambitious policy goals to protect biodiversity and maintain the resilience, functions, and services that ecosystems provide (Balmford et al., 2010;Steffen et al., 2015;IPBES, 2019;Hermoso et al., 2022).
To overcome this shortcoming, we tested an indicator-based approach for assessing ecological resilience by estimating, comparing, and visualizing a suite of potential resilience indicators on the basis of readily available monitoring data from the North Sea. The North Sea is a heterogeneous marine ecosystem with large regional variations in environmental conditions and human impacts, including depth, temperature, primary production, and fishing effort (Reiss et al., 2010;Couce et al., 2020). These conditions profoundly affect individual species distribution and abundance, as well as the overall composition and structure of fish communities. The primary spatial pattern identified among the set of indicators is a marked north-south gradient. In general, we found that the northern central part of the North Sea is characterized by a higher number of functional groups but generally lower redundancy and evenness within functional groups. This spatial structuring generally conforms with findings from previous studies on fish community traits and diversity (Reiss et al., 2010;Dencker et al., 2017;Pecuchet et al., 2017;. These studies demonstrate that communities in the northern part of the area are primarily composed of large, latematuring, long-living, and slow-growing species that feed at high trophic levels. On the contrary, communities in the southern parts of the North Sea are dominated by smaller, fastergrowing species feeding on lower trophic levels (Dencker et al Figures S6 and S7 for corresponding plots based on the lower and higher levels of grouping using 10 and 50% cutoffs.) The metrics used in this study represent a range of functional diversity indices that can uncover the complex natural structures of communities. Using multiple metrics allows for estimating and visualizing key components of resilience based on long-term survey data. However, it is important to explore different types of indicators and metrics, as they can be context and system dependent. In this study, we included multiple metrics to explore and compare the consistency of the resilience indicators. For functional group evenness, we observed a strong consistency between the two selected metrics (EvennessSh and EvennessPE). Consequently, future applications of our indicator-based approach could consider using only one of the evenness metrics as a resilience indicator. For further analysis, we continued with EvennessSh. However, in terms of functional response diversity, the individual metrics do not show the same robustness, as there was quite some spatial heterogeneity and differences between the two metrics (response diversitysyn and response diversitypoe). Al-though metrics of response diversity have a great potential for representing key components of resilience and stability (Loreau and Mazancourt, 2008;Anderson et al., 2013;Lindegren et al., 2016), they require persistent and high-quality temporal data on species abundances. Consequently, they are sensitive to sporadic, rare, and less abundant species whose true abundance fluctuations may not be reliably represented by the sampling method. To that end, caution should be exercised when using these indicators and metrics, conditioned on the frequency and resolution of the underlying sampling. Other relevant measures of functional diversity could also be explored in future applications, such as functional attribute diversity (Walker et al., 1999), functional divergence (Villéger et al., 2008), and Rao's quadratic entropy (Rao, 1982;Botta-Dukát, 2005). These additional measures would fall within the broader four categories of indicators in our framework and could be readily included in a resilience assessment. Furthermore, we recommend investigating temporal trends and changes in individual resilience indicators, at least for metrics not requiring long-term time series for their estimation.
In addition to the observed coherence or differences between individual metrics, our study shows similarities and differences between the selected set of indicators. As an illustrative example, functional redundancy and evenness both load positively on PC1 of our PCA analysis, while the indicators of functional group richness and response diversity are highly correlated with PC2. These similarities and differences are, in turn, evident from the individual spatial patterns of each indicator. Since the PC axes are orthogonal to one another, this shows that the set of indicators and metrics is representing different and complementary aspects of resilience. Hence, accounting for a range of indicators is important to cover different aspects potentially contributing to resilience. Furthermore, our results imply that high resilience cannot be achieved by maximizing all individual aspects of resilience, simply because there seem to be inherent trade-offs between these components. Notably, high functional group richness, as observed in the northern North Sea, does not necessarily lead to higher redundancy, nor evenness among species sharing similar traits. This is likely due to certain functional groups being occupied by only a single or handful of species (Mouillot et al., 2014;Lindegren et al., 2018), potentially with very different relative abundances (leading to low evenness). To that end, we stress the need for future work on functional ecology to explore the dependence and potential trade-offs between richness, redundancy, evenness, and response diversity affecting community properties. Since such trade-offs can be system dependent, we therefore recommend future empirical studies to investigate the influence of resilience indicators on overall system resilience, stability, and functioning.
To facilitate such future investigations and insights, it is important to formulate a set of expectations and hypotheses regarding the links and causal relationships between resilience indicators and the overall properties of ecosystems. Below, we briefly contextualize the spatial patterns in resilience indicators by elaborating on their potential roles and mechanisms supporting ecosystem functioning and stability at large, with reference to relevant literature. First of all, the choice of indicators and their associated metrics included in this study build on already established resilience theory, which indicates that increased biodiversity ensures ecosystem stability, functions, and services (Buchmann, 2006;Cardinale, 2006;Cardinale et al., 2012). Hence, the selected metrics constitute a suite of biodiversity measures shown to be important determinants of ecological resilience (Hillebrand and Matthiessen, 2009;Angeler and Allen, 2016;Sundstrom et al., 2018). High levels of richness have been demonstrated to increase ecosystem functioning and stability over time (Tilman et al., 1996;Hillebrand and Matthiessen, 2009). Moreover, ecosystem functions, such as flux of matter, energy, and primary production, depend on functional attributes of species present in the ecosystem (Cadotte et al., 2017). Functional composition and diversity are therefore considered more critical for ecosystem processes and resilience within a system than the number of species itself (Gagic et al., 2015). Redundancy occurs when multiple species share similar traits and functional roles. This provides critical insurance and functional reinforcement (Walker, 1995;Walker et al., 2004). Consequently, as species richness increases, there will be an overlap in the functional traits among species within a given community. This niche overlap secures ecosystem functioning and stability over time because if one species is lost, a similar species may serve to retain the overall niche and functionality (Sundstrom et al., 2018). In addition to the richness and redundancy of species within functional groups, the relative abundances of species may regulate ecosystem functions and processes (Hillebrand and Matthiessen, 2009;Maureaud et al., 2019). Notably, it has been suggested that high evenness allows communities to quickly adapt to new environmental conditions and maintain high production over time if exposed to changes (Hillebrand and Matthiessen, 2009). Lastly, ecological communities can stabilize ecosystem processes if species vary in their responses to environmental fluctuations such that an increased abundance of one species can compensate for the decreased abundance of another (Yachi and Loreau, 1999;Gonzales and Loreau, 2009). Such response diversity, often characterized by asynchrony in species fluctuations within functional groups (Lindegren et al., 2016), can support ecosystem resilience in the face of anthropogenic pressures and environmental uncertainties, as it refers to the range of reactions to environmental change among species that contribute to the same ecosystem function (Elmqvist et al., 2003). To distinguish regional differences in the potential overall resilience, we clustered out four areas based on the indicators of resilience. This information and their potential links to ecosystem functions and stability may shed light on geographical differences in the overall resilience and sensitivity of these communities to natural and human pressures. This can provide important information to regional ecosystem assessments and advice. For instance, the resilience indicators and maps may be included in the ICES Integrated Ecosystem Assessments (notably WGINOSE) as well as the associated ecosystem overview documents supporting evidencebased ecosystem-based management. The indicators and maps may also be taken up within existing regional seas conventions, such as OSPAR and HELCOM, where they may complement their current list of indicators used in assessing biodiversity, or good environmental status, broadly. Furthermore, the detailed maps of resilience indicators, especially if estimated and provided at an even finer spatial resolution, may be more directly applicable in vulnerability assessments and marine spatial planning efforts (Gladstone-Gallagher et al., 2019;Harris et al., 2022), notably, in order to prioritize the distribution of marine protected areas to safeguard and restore biodiversity, as called for by the EU Green Deal and the associated EU Biodiversity Strategy 2030. To that end, accounting for such structuring may secure the integrity of the system, which provides ecological resilience and ensures a dynamic community that can adapt to environmental changes and threats (Holling, 1973). To do so, it is important to provide tools that are capable of representing ecological resilience and its changes in space and time. The resilience indicators applied in this study are generic and relate to fundamental properties of an ecosystem. The indicators can be estimated from already existing data, which makes it possible to measure and monitor spatio-temporal changes in largescale ecosystems. This is getting increasingly easier thanks to the availability of large datasets from long-term monitoring, but further accessibility is still essential (Maureaud et al., 2021). Furthermore, up-to-date and accessible trait data on species are necessary to cover the true diversity of species. This would allow for a diverse selection of species traits covering the ecological niche of the species and their function.
Due to the difficulty of directly measuring resilience of ecosystems per se, our study shows that ecological resilience can instead be described by a suite of complementary functional diversity indicators and metrics resembling the underlying attributes promoting resilience. Hence, we argue that this indicator-based framework has great potential to represent the resilience of ecosystems to change and may serve to inform management and conservation actions in marine, terrestrial, and freshwater systems.