Evidence of climate‐driven regime shifts in the Aegean Sea’s demersal resources: A study spanning six decades

Abstract Climate change (CC) can alter the configuration of marine ecosystems; however, ecosystem response and resilience to change are usually case‐specific. The effect of CC on the demersal resources of the Aegean Sea (east Mediterranean Sea) was investigated during the past six decades applying a combination of multivariate analysis, non‐additive modeling and the Integrated Resilience Assessment (IRA) framework. We focused on the study of: (i) the biological “system” complex, using proxies of biomass (landings per unit of capacity) for 12 demersal taxa, and (ii) the environmental “stressor” complex, described by 12 abiotic variables. Pronounced changes have occurred in both the environmental and biological system over the studied period. The majority of the environmental stressors exhibited strikingly increasing trends (temperature, salinity, primary production indices) with values started exceeding the global historical means during late 1980s‐early 1990s. It is suggested that the biological system exhibited a discontinuous response to CC, with two apparently climate‐induced regime shifts occurring in the past 25 years. There is evidence for two‐fold bifurcations and four tipping points in the system, forming a folded stability landscape with three basins of attraction. The shape of the stability landscape for the Aegean Sea's biological system suggests that while the initial state (1966–1991) was rather resilient to CC, absorbing two environmental step‐changes, this was not the case for the two subsequent ones (intermediate: 1992–2002; recent: 2003–2016). Given the current trajectory of environmental change, it is highly unlikely that the biological system will ever return to its pre‐1990s state, as it is entering areas of unprecedented climatic conditions and there is some evidence that the system may be even shifting toward a new state. Our approach and findings may be relevant to other marine areas of the Mediterranean and beyond, undergoing climate‐driven regime shifts, and can assist to their adaptive management.


| INTRODUC TI ON
Over 90% of the energy gained by the Earth as a result of human activities accumulates in the ocean, and most of it stays in the upper 700 m of the water column (Rhein et al., 2013). Hence, the oceans "feel" most of the impacts of climate change (CC). The impacts are not geographically homogeneous, with the North Atlantic and Mediterranean Sea storing a disproportionally large stock of anthropogenic carbon compared to other ocean regions (Khatiwala et al., 2013). Climate variability and change affect marine species in a multitude of ways from eliciting physiological responses and changing productivity and survival, to altering habitat and resource availability that lead to range expansions/contractions and shifts of geographical distributions (Koenigstein et al., 2016).
Apart from direct effects, a change of abiotic conditions may affect species indirectly through their interactions with other species.
Due to the interconnectedness of biological systems, CC may trigger a complete reorganization of an ecosystem, known as "regime shift" (Scheffer, 2009;Scheffer & Carpenter, 2003). Regime shifts are described as abrupt changes of the ecosystem state (e.g., characterized by different dominance patterns) and can be triggered even by smooth, continuous climatic shifts, particularly in ecosystems with low ecological resilience, for example, as a result of anthropogenic pressures (fishing, pollution, habitat change). Typically, after a regime shift, an ecosystem state is self-stabilizing through positive feedback loops and the restoration of the environment to previous conditions does not necessarily cause a recovery of the ecosystem to its previous state (hysteresis) (Scheffer, 2009). In this context, alternate regimes refer to dynamic states that can be viewed as basins of attraction around the system response curves (attractors), forming a folded stability landscape on the system-stressor surface (Barnosky et al., 2012;Scheffer, 2009)-see Figure S0. Systems within a folded stability landscape may switch to a new state either through a "horizontal movement" across a tipping point associated with a change in stressors, or through a "vertical movement" associated with a small change in the system state while the system lies close to a tipping point (Scheffer, 2009).
Two major climatic shifts have been documented in the North hemisphere during the 20th century: in the mid-1970s and the late 1980s (Möllmann & Diekmann, 2012 and references therein). Both were related to a change of the North Atlantic Oscillation (NAO) Index, which was in turn associated with an increase in sea temperature. Concurrent with these climatic shifts, several north European marine ecosystems (North Atlantic, North and Baltic Seas, Scotian shelf, and Black Sea) were reported to have undergone ecosystem regime shifts (Möllmann & Diekmann, 2012). All these studies suggest that besides fishing, CC is a significant factor shaping the fluctuations of marine resources (Barange et al., 2018).
In the Mediterranean Sea, a regime shift in oceanographic variables and planktonic assemblages has been documented in the western part of the basin in the late 1980s (Conversi et al., 2010).
Additionally, regime shifts in fisheries resources have been reported across the Mediterranean basin over the past 30 years (Tzanatos et al., 2014), with such shifts being the result of ecosystem reorganizations as a discontinuous response to sea warming (Vasilakopoulos et al., 2017). Furthermore, the exploitation rate of Mediterranean stocks over the past 20 years has been steadily increasing, selectivity (proportional exploitation of juveniles) has been deteriorating, and stocks have been shrinking (Vasilakopoulos et al., 2014). It is likely that Mediterranean fisheries vulnerability to CC is higher, given overfishing, higher exposure to warming, arrival of non-indigenous species, and an overall lower adaptive capacity (Hidalgo et al., 2018).
Herein, we combine multivariate and nonadditive modeling tools to identify plausible regime shifts in the demersal community of the Aegean Sea, implementing the Integrated Resilience Assessment (IRA) framework (Vasilakopoulos & Marshall, 2015;Vasilakopoulos et al., 2017). We apply this framework to an empirical biological data-series (fisheries-dependent relative biomass) during 1966-2017, and investigate the effect of the concurrent development of a series of environmental variables. The analysis of this exceptionally long time-series allows us to unveil previously unknown climatedriven ecosystem dynamics in the Aegean Sea and may assist in future adaptive management of the local marine biological resources.

| Study area-fisheries
The study area covers the Aegean Sea (Figure 1), located in the Eastern Mediterranean Sea. Generally, the region has a complex geomorphology, which directly or indirectly influences the fishing activity; the Aegean can be divided in two distinctive ecoregions: the south Aegean is an oligotrophic area where evaporation exceeds freshwater income, while the north Aegean is influenced by the cold, brackish, nutrient-rich Black Sea waters (Lykousis et al., 2002;Zervakis et al., 2004). As a result, more than 50% of the Greek otter bottom trawl fleet targeting demersal species, operates in the north Aegean Sea, producing >57% of the total national demersal landings; north Aegean bottom trawl landings exceed those in the south Aegean by a factor of more than 5 (HELSTAT, 1964(HELSTAT, -2017. Approximately 15,000 commercial vessels are involved in the Greek demersal fisheries (using primarily bottom trawls, static nets and longlines), 12,000 of which are exerting their effort in the study area (Aegean Sea) (HELSTAT, 1964(HELSTAT, -2017. The Greek demersal fisheries are multispecies in nature targeting a highly diversified mix of fish, cephalopods, and crustaceans. Management is conducted through input controls (spatio-temporal effort restrictions) and various technical measures (Kapantagakis, 2007).

| Environmental
A large number of historical physicochemical variables for the Aegean Sea during the period 1966-2017 were compiled (Table S1).
Data were made available through the Horizon2020 CERES project (http://ceres proje ct.eu/). They were derived from the POLCOMS-ERSEM coupled model framework, producing daily, monthly, and annual outputs of physical, chemical, and biological variables at a 1/10 of a degree spatial resolution. POLCOMS (Proudman Oceanographic Laboratory Coastal Ocean Modelling System, Holt et al., 2001) is a 3-D community physical model with the ability to run in regions that include both the deep ocean and the continental shelf. The European Regional Sea Ecosystem Model (ERSEM; Baretta et al., 1995) is one of the most complex lower trophic-level marine ecosystem models currently in use, including bacteria, four phytoplankton and three zooplankton functional groups, a fully resolved diurnal cycle, variable carbon to chlorophyll ratios and independent nutrient pools for carbon, nitrogen, phosphorous, and silicate. The POLCOMS-ERSEM modeling system is well established in the NE Atlantic and has been applied successfully in the Mediterranean Holt et al., 2001Holt et al., , 2012 To shortlist the environmental variables of interest to this study, we focused on their relevance for the demersal resources ("system"). From the full set of 36 environmental variables, some refer to the sea bottom (e.g., ben_DOC-benthic dissolved carbon), others are by definition 2D (e.g., optical depth), while most of them were available at different depth zones (3D variables, e.g., gross primary production). From a series of experimental surveys in the Aegean Sea (Katsanevakis et al., 2009), most species/taxa analyzed herein (see Demersal resources below) exhibited higher abundance at depths of <50 m. Hake differentiated from this general pattern, showing high abundance both in waters beyond the continental shelf (Katsanevakis et al., 2009;Labropoulou et al., 2008;Sion et al., 2020) as well as in shallower waters (PANDORA, 2019).

F I G U R E 1 Bathymetric map of the study area
Hake abundance is linked to the occurrence of chlorophyll-a fronts and high densities are to be expected after local phytoplankton blooms (Cartes et al., 2004). Alemany et al. (2014) andDruon et al. (2015) showed that surface oceanic features are strong drivers even for demersal resources, such as hake. Sion et al. (2020) found that Mediterranean hake biomass indices exhibited similar trends for sea bottom temperature and sea surface temperature.
Sea surface temperature is typically used to express the preferred temperature of demersal species (e.g., Cheung et al., 2013 and subsequent literature on the Mean Temperature of the Catch-MTC) and to capture the effect of warming upon marine ecosystems, including demersal resources, within Integrated Ecosystem Assessments (IEA; e.g., Möllmann & Diekmann, 2012 and references therein). Moreover, employing correlation analysis on the available variables by depth, confirmed that the concern over the choice of depth can be relaxed, since their values at different depths were highly correlated ( Figure S2). Consequently, for 3D variables it was decided to use yearly averages at 5 m depth over the whole Aegean Sea area. To avoid duplicating information, correlation testing (Szekely et al., 2007) was used to remove highly correlated variables, ending up with a subset of 12 variables ( Figure S2; Figure S1): T, S, ssh, grossPP, pH, O, ben_DOC, MLD, Chl1, Chl2, Chl3, P3C (full description in Table 1). The use of HP or kW as a surrogate for effort would have led to unreliable estimates, due to the renowned underdeclaration of engine power (Damalas et al., 2014;DG MARE, 2019); as a result, Gross Tonnage (GT) was used as the most consistent capacity measure of the fleet targeting demersal stocks. GΤ has been successfully applied in an analogous long-term study in the Mediterranean Sea (Fortibuoni et al., 2017).

| Demersal resources
Landings of demersal species in the Aegean Sea account for ~25% of total Greek landings (HELSTAT, 1964(HELSTAT, -2017. In this study, 12 demersal taxa were considered, based on the availability of uninterrupted time series, and were expressed in Landings Per Unit of Capacity (LPUC), defined as the sum of annual landings of each species from the demersal fleets (trawlers and smallscale fishery) divided by the sum of GT of these fleets (kg/GT).
To confirm that this GT-based LPUC was a good proxy of biomass, we investigated the correlation between LPUC and LPUE (effort in GT x DAS, kW x DAS) for the most recent years (e.g.,

>2002
) that detailed effort data were available through the EU DCF. A statistically significant positive correlation was confirmed (p < .05- Figure S3). The period of data used spans more than half a century (from 1966 to 2017; Figure 2).

| Statistical analyses
The IRA framework was used to study the multivariate response of the demersal resources ("system") to environmental change ("stressors") and carry out a resilience assessment (Vasilakopoulos & Marshall, 2015;Vasilakopoulos et al., 2017). The IRA is an extension of the IEA framework (Levin & Möllmann, 2015;Möllmann & Diekmann, 2012;Möllmann et al., 2009). On top of the IEA routines used to study the multivariate development of complex natural systems, the IRA applies non-additive modeling on the system-stressor relationship and quantifies ecological resilience to construct stability landscapes.
As a first step, we reduced the complexity of the multivariate datasets by means of a multivariate analysis (principal component analyses, PCA) and implemented a suite of IEA methods to identify

Variable
Units Description The analysis was applied to the biological system and environmental stressor datasets separately, to find holistic indicators of both the system's and stressors' state and identify possible alternate regimes.
Subsequently, to investigate the system-stressor relationship, the system holistic indicator was regressed against the stressors' indicator using additive (continuous) statistical models (GAMs). The fit of the GAM models was then compared to the respective fit of non-additive (discontinuous) statistical models (TGAMs) to assess the type (continuous/discontinuous) of the system response and estimate the potential alternate attractors and threshold years of the system. Finally, a resilience assessment was carried out and a stability landscape was constructed. The analysis was carried out in R (R Core team, 2019) using packages energy ( can be found in Rao (1964) and Legendre and Legendre (1998).
The temporal development of each variable in the data set was according to their loadings on the PC1 axis, raw values of each variable were categorized into quintiles, and each quintile was given a specific color. For each variable, values of the lowest and highest quintile were indicated in green and red, respectively, with another three levels of intermediate colors.

| Chronological clustering
Chronological Clustering is a variation of cluster analysis, designed to run on ordered/successive samples like time-series (Andersen et al., 2009;Legendre et al., 1985). The method calculates a similarity matrix combining samples into clusters, based on a permutation test and predefined significance and connectedness levels (Legendre et al., 1985). In our case, chronological clustering was applied using the "coniss" method (Constrained Incremental Sums of Squares- Grimm, 1987) on the 12 demersal species/stocks LPUC time-series, as well as on the 12 environmental variables.

| Breakpoint structural analysis
Breakpoint analysis (Zeileis et al., 2002) allows identifying statistically significant changes in the mean level of subsets of a time series.
A time series is randomly split in two or more subsets (''data windows'') and the mean levels are compared through a modified F test (''structural change'' or sctest- Zeileis et al., 2010). The procedure is repeated iteratively until all significant breakpoints (if any) are identified. The Bayesian Information Criterion is used as an objective criterion to determine the number of breakpoints and their associated dates with the corresponding 95% confidence intervals. In this study, breakpoint structural analysis was applied to the first two principal components of PCAsys and PCAstr to assess the year(s) of statistically significant change. All analyses were performed in R library strucchange (Zeileis et al., 2002).

| STARS
The Sequential t-Test Analysis of Regime Shifts (STARS) is based on a sequential Student's t-test that signals the possibility of a regime shift (Rodionov, 2004); the identified year of regime shift, is indicated by the calculated probability level. The method consists of calculating a Regime Shift Index (RSI), to accept or reject the hypothesis of a regime shift at each observation. RSI represents a cumulative sum of normalized anomalies relative to a critical value (Rodionov, 2006). The time scale of shifts to be considered is controlled primarily by the cut-off length, which is similar to the cut-off point in lowpass filtering, and determines the minimum length of the regimes for which the magnitude of the shifts remains intact (Rodionov & Overland, 2005). A longer cut-off length hence identifies the strongest signal (as opposed to many smaller events). Herein, STARS was applied to detect significant shifts in the time-series of the first two principal components of PCAsys and PCAstr during 1966-2017 using a cut-off length of 15 years. A tool for the detection of a regime shift based on STARS is available on the website http://www. berin gclim ate.noaa.gov/regimes.

| TGAMs
TGAMs are special cases of GAMs (Generalized Additive Models- Hastie, 2017). GAMs 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) signifying that a relationship is best described by two different functions, one before and one after the threshold (Casini et al., 2009;Ciannelli et al., 2004). These different functions are visualized as separate response curves (TGAM branches).
Threshold GAMs were used to regress environmental variables against biological variables (LPUC) and identify different response "regimes." To compare the fit of continuous (GAMs) and discontinuous models (TGAMs) on the relationships between the system indicator PC1sys and environmental stressors PC1str, at different time lags, the "genuine" cross-validatory squared prediction error (gCV) was computed (Ciannelli et al., 2004). The gCV accounts for the es- It should be noted that while chronological clustering, breakpoint analysis, and STARS detect shifts in the relationships between each of the holistic system and/or stressor variables and time, TGAMs detect shifts in the relationship between the holistic system and a stressor variable. Hence, it is possible to have shift(s) in the system and/or stressors identified by chronological clustering/breakpoint analysis/STARS, but no shift in the system-stressor relationship (in which case a GAM would provide a better fit than a TGAM), and vice versa. For example, if the system shifts at the same time when the stressors shift, then the system-stressor relationship could be linear. However, a shift in the system that coincides with a shift in the system-stressor relationship, suggests a system shift, that is, a discontinuous response to changing stressors.

| Resilience assessment
Following the IRA framework (Vasilakopoulos et al., 2017), the branches of the optimal TGAM model(s) were considered to represent the alternate attractors of the system, forming foldbifurcations. The plausibility of the transition trajectories (i.e., shifts happening at local extrema of stressors followed by hysteresis) was then confirmed (Tsimara et al., 2021). Resilience is expected to increase with increasing distance of a system state to the tipping point ("horizontal component" of resilience-hComp; Figure   S0). In contrast, resilience is expected to decrease as the distance of the system state to its attractor increases ("vertical component" of resilience-vComp; Figure S0). The annual resilience estimate (Res y ) was thus the sum of the horizontal distance

| Environmental "stressors"
Visual inspection of the anomalies (
We cannot be conclusive regarding the last shift (2016)

| System-stressors relationship
TGAMs offered superior fits to GAMs when modeling the regression of the biological "system" indicator (PC1sys) against the environmental "stressors" (PC1str) ( Table 2). In all different time lags investigated, 1991 was found to be the threshold year. The 0-lag model was found to provide the best fit amongst the tested TGAMs, but the difference in gCV between models with different lags was quite small (Figure 9; Table 2 The outcomes of all the methods employed herein are summarized in Table 3. There is a striking convergence of all approaches to almost the same years suggested as threshold years; an indication of a clear signal in the data that two regime shifts occurred in the system as a discontinuous response to changing environmental stressors.

| Folded stability landscape
The resilience assessment allowed the construction of a folded stability landscape of the biological system that exhibited three basins of attraction (regimes) (Figure 10). It should be noted that the area between tipping points (F1 and F2; F3 and F4) hosts barely any observations, as it is an area of repulsion; this is a common feature of stability landscapes. Hence, the linear interpolation of resilience between the years of the filled contour often creates "corridors" of F I G U R E 3 Anomaly of environmental variables expressed as difference from the overall mean divided by the standard deviation in the Aegean Sea during 1966-2017 seemingly higher resilience between years with high resilience belonging to different states (Vasilakopoulos et al., 2017). To correct for this, we added two theoretical "years" (not plotted) between F1 and F2, and between F3 and F4, whose resilience was calculated based on their distance from the attractors and tipping points, exactly like for all other years.
The folded stability landscape suggested that the system retained its original state for 25 years despite concurrent large changes in the environment. The system even absorbed two environmental step-changes without changing its state, owing to its high resilience. However, the system's resilience was gradually eroding in response to environmental change, and in 1990-1992, a small perturbation,

| D ISCUSS I ON
A combination of multivariate analyses, non-additive modeling, and a resilience assessment, that is, an IRA, was employed to investigate a plausible effect of CC on the demersal resources of the Aegean Sea (east Mediterranean Sea) over the past six decades. This analysis revealed that pronounced changes have occurred in both the environment and the biological system over the studied time period. It is suggested that the biological system has exhibited a discontinuous response to CC, with two climate-induced regime shifts occurring in the past 25 years. There is evidence for twofold-bifurcations and four tipping points in the system, forming a folded stability landscape with three basins of attraction (Scheffer, 2009).
The majority of the environmental stressors investigated exhib-  Figure 10) compared to the initial state. It could be that environmental change, up to 1991, represented "familiar" conditions that had occurred during the evolutionary history of the biological system; hence, the system had built some resilience to them.
However, after 1991, the environment might have moved to areas never experienced before; hence, the system struggled to stabilize within its new states. The likely onset of a new state in 2016/17 corroborates this hypothesis. The biological system seems to be increasingly prone to discontinuous regime shifts in response to CC over the past 25 years as it is entering areas of unprecedented climatic conditions. This is an important finding that might be also relevant to other marine areas moving into novel climatic regimes as a result of the global CC.
At a global scale, changes in marine fisheries catch have been documented to be related to changes in ocean features (e.g., temperature), with an increasing dominance of catches of warmer waters species at higher latitudes and a corresponding decrease in the proportion of catches of subtropical species in the tropics (Cheung et al., 2013). Marine biota responds to ocean warming through changes in distribution, abundance, phenology, and body size, leading to alteration of community structure and trophic interactions ultimately affecting fisheries catches (Cheung et al., 2010). In general, species and marine communities respond differently to change, depending on the organisms' capacity to adapt to novel regimes.
Lower trophic level organisms are more able to integrate seasonal thermal changes, whereas, higher trophic level species depend on both the availability of lower trophic level organisms and environmental conditions (Friedland et al., 2019).
The Mediterranean Sea is a semienclosed basin, sensitive to bottom-up processes (Piroddi et al., 2017), where the effects of CC are likely to be more pronounced and become apparent sooner than in other, more open ocean regions (Giorgi, 2006). The Mediterranean basin is projected to be a hotspot of CC in terms of both warming and an increase in extreme events (Calvo et al., 2011;de Madron et al., 2011;MEDCLIVAR, 2004). Since the late 1990's, the abrupt increase in temperature has facilitated the establishment and northward migration of warm-water species from adjacent seas at an unexpectedly rapid rate (Raitsos et al., 2010).
Within the Mediterranean, the Aegean Sea in particular, encompasses certain peculiarities that may categorize it among the most vulnerable regions: an extended length of coastline (12,000 kmalmost 1/4 of the whole Mediterranean coastline), complex bathymetry, and numerous islands (more than 3000 islands or islets), hosting various intensive anthropogenic activities such as industry, aquaculture, maritime shipping (over 35,000 vessels per year traversing through its waters), and tourism (e.g., 800,000 yachts sailing in the Mediterranean, Aegean Sea being one of the most popular sailing destinations), (Sakellariou et al., 2005;Sofianos et al., 2002).
However, its most striking characteristic is the largest fishing fleet in Europe (~12,000 vessels, out of a total of 15,000 in Greece and 80,000 in the EU) 1 , mostly comprising of small-scale fishery vessels, operating in a traditional manner and being multi-specific in nature (HELSTAT, 1964(HELSTAT, -2017. Adding to this vulnerability, the region has been probably extensively exploited by humans earlier than any other marine region of the world (NHRF, 2010). The most recent warning sign comes from the acceleration in the annual mean rate of warm and tropical alien species introduced (a 150% increase of species entry in the last two decades), with the speed of alien species spreading being faster than environmental change (EASIN, 2020; Raitsos et al., 2010). In the Mediterranean, Aegean Sea is now considered among the regions being hot-spots for alien species (EASIN, 2020;Zenetos et al., 2011).
The potential vulnerability of the fishing sector in the Aegean Sea to the effects of CC needs to be better understood and considered to inform the development of effective management and adaptation policies. Free et al. (2019) highlight the importance of  (Lloret et al., 2015). Future projections for the Mediterranean predict that the western basin is expected to become more oligotrophic with a surface density decrease, while surface production is expected to increase in the eastern basin linked to a density increase (Macias et al., 2015). Regional changes in fish abundance and their distribution will alter species richness, with an expected increase in overall richness by the mid-21st century in the Eastern Mediterranean, and a decrease in the western region (Albouy et al., 2013). This means that, while the biological transitions in the Aegean Sea detected here are of profound importance for local fisheries management, which should be tailored to the most recent ecosystem state, we should also prepare for a future where abrupt ecosystem shifts to novel states will become the norm.
Consequently, managers and fishers need to focus on increasing the adaptive capacity of the fishing sector to protect its socio-economic resilience in the context a decreasing ecological resilience. This task could encompass diversification of the fishing activities, adjustment of the timing and location of spatio-temporal measures, as well as new market regulations.
While our study has focused on the effect of CC, other potential stressors of the biological system may be at play, with the most important one being fishing. The absence of key fishing data for the most part of the time series (i.e., before 2000; e.g., actual fishing effort, stock status) hindered a deeper understanding of potential fishing effects. Nevertheless, a LPUC index was used as a proxy for biomass with vessel capacity (in GT) as a surrogate of fishing effort.
Based on these limitations, certain assumptions had to be made: (i) fishery performance can approximate, reasonably, the relative biomass at sea (Marr, 1951), and (ii) commercial fishery statistics is a fishery performance index (Ricker, 1975). Obviously, other fishery effects (e.g., technological creep) and other external drivers may exist interacting in ways not fully identified herein.
The Greek fleet has undergone some dramatic changes during the study period; Greece joining the EU in the early 1980s, gave access to funds allowing for investing in and expanding the fleet.
However, the EU Common Fisheries Policy reforms in 1992 and 2002 (EC, 2002a(EC, , 2002bEEC, 1992) requesting adjustment of fishing capacity of the fleets, led to subsidies for mass decommissioning and scrapping of vessels. As a result, the Greek fleet contracted by more than 30% in just 20 years, from more than 22,000 vessels in 1995 down to almost 15,000 vessels in 2016 ( Figure S5). Yet, the effect of fishing has been suggested as a driver of stock status only during the past two decades (e.g., Vasilakopoulos et al., 2014 (1991) and stressor systems provides strong evidence that CC has been a major driver and that the human factor alone could not lead to such a conspicuous event. The main shift identified herein (1992) occurred during a period that fishing capacity (and presumably fishing intensity) remained stable (1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996); Figure S5). This strengthens the argument that the key driver behind the observed system dynamics is CC rather than fishing. Nevertheless, it could be argued that this sensitivity of the system to CC was related to some extent to the overfished status of the fisheries resources (Hsieh et al., 2006).
Further understanding of the interactions and synergies between CC and fishing, and the ways they affect different trophic levels and the ecosystem as a whole is needed.

Method
Estimated threshold year(s) Environmental "stressors" Biological "system" Chronological clustering 1987 (secondary cluster split at 1969) 1990 (secondary cluster split at 2005) Breakpoint analysis 1982, 19881990STARS 1975, 1987(high overlap) 1991 TGAMs (PCA1sys ~ PCA1str) 1991, 2002 Note: First three methods identify shifts in the relationship of the composite environmental "stressors" and biological "system" and against time, while TGAM identifies shifts in the relationship between the biological "system" and environmental "stressors." F I G U R E 1 0 Empirical folded stability landscape with transition pathways of the Aegean demersal landings assemblage ("LPUC profile") in response to environmental change during the period 1966-2017. The folded stability landscape exhibits two basins of attraction and four tipping points (F1, F2, F3, F4). Color scale indicates the relative resilience ("altitude"), continuous lines indicate the attractors, dotted black lines indicate the possible extensions of the attractors, and gray dashed lines indicate the borders of the basins of attraction. The prefix 19-and 20-of the years is omitted. Year 2017 has not been assigned to any basin of attraction as it differed substantially from its preceding years dominated by cephalopods and crustaceans, in an environment characterized by elevated water temperatures, higher salinities and more productive and acidic waters. Furthermore, the system needs to be closely monitored as it moves toward an area of unprecedented environmental conditions; a terra incognita, which could entail further discontinuous responses to CC. Resource Assessments).

CO N FLI C T O F I NTE R E S T
None declared. Funding acquisition (equal); Writing-original draft (supporting);

AUTH O R CO NTR I B UTI
Writing-review & editing (supporting).

DATA AVA I L A B I L I T Y S TAT E M E N T
Data on historical physicochemical variables for the Aegean Sea during the period 1966-2017 were made available through the Horizon2020 CERES project (http://ceres proje ct.eu/). Marine fisheries' landings and fleet capacity data for the Greek fleet have been extracted from the yearly bulletins of the Hellenic Statistical Authority (HELSTAT, 1964(HELSTAT, -2017