# Detecting Regime Shifts in the Portuguese Continental Shelf Ecosystem Within the Last Three Decades

Dorota Szalaj, Alexandra Silva, Pedro Ré, Henrique Cabral

► To cite this version:

Dorota Szalaj, Alexandra Silva, Pedro Ré, Henrique Cabral. Detecting Regime Shifts in the Portuguese Continental Shelf Ecosystem Within the Last Three Decades. Frontiers in Marine Science, 2021, 8, pp.24. 10.3389/fmars.2021.629130. hal-03290115

HAL Id: hal-03290115

https://hal.inrae.fr/hal-03290115

Submitted on 19 Jul 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for
updates

# Detecting Regime Shifts in the Portuguese Continental Shelf Ecosystem Within the Last Three Decades

Dorota Szalaj<fnref n="1" />,<fnref n="2" />,³*, Alexandra Silva<fnref n="1" />,<fnref n="2" />, Pedro Ré<fnref n="1" /> and Henrique Cabral<fnref n="1" />,⁴

Marine ecosystems are affected by diverse pressures and consequently may undergo significant changes that can be interpreted as regime shifts. In this study we used integrated trend analysis (ITA) that combines multivariate statistics and methodologies to identify abrupt changes in time-series, in order to test a hypothesis about the occurrence of regime shifts in the Portuguese continental shelf ecosystem (PCSE). We used two types of data describing ecosystem drivers (fishing mortality and environmental/climatic indices) and ecosystem state (observed and modelled biomass and ecosystem indices). Modelled biomass and ecosystem indices were outputs of Ecopath with Ecosim temporal model parametrised for PSCE between 1986 and 2017. The analyses indicated that the regime shifts in the PCSE have occurred during three periods in the last three decades: “early regime” until the mid-1990s, followed by “transition regime” in-between and “late regime” since the mid-2010s. The detected regime shifts are characterised by changes in the pelagic community that became more dominant when compared to the demersal community and shifted from sardine, the main fishing resource, abundant in the “early regime”, to other less valuable pelagic fishes such as chub mackerel that dominated the “late regime”. The “early regime” was characterised by high catch, a larger proportion of demersal species, and higher diversity while, the “late regime” was represented by lower catch, an increase in higher trophic level (TL) predatory fish and lower diversity. Moreover, the “late regime” showed lower resilience and reduced maturity when compared to the “early regime”. Changes described in the ecosystem were probably related to (1) the shift in the north Atlantic environmental conditions that affected small pelagic fish (SPF) and lower TLs groups, (2) reduction in fishing pressure, and (3) internal triggers, related to the indirect trophic interactions that might have benefited higher TL fish and impacted the pelagic community. In the context of PCSE management, this study highlighted a need to consider the possibility of regime shifts in the management process. For example, regime specific harvest rates and environmental reference points should be considered when an indication of abrupt change in the ecosystem exists.

**Keywords:** Portuguese continental shelf ecosystem, regime shift, ecosystem modelling, integrated trend analysis, environmental forcing, sardine decline

**OPEN ACCESS**

**Edited by:**

Aschwin Hillebrand Engelen,
University of Algarve, Portugal

**Reviewed by:**

Francisco Leitão,
University of Algarve, Portugal
Maciej Tomasz Tomczak,
Stockholm University, Sweden

**\*Correspondence:**

Dorota Szalaj
dszelaj@fc.ul.pt

**Specialty section:**

This article was submitted to
Marine Ecosystem Ecology,
a section of the journal
Frontiers in Marine Science

**Received:** 13 November 2020

**Accepted:** 05 March 2021

**Published:** 26 March 2021

**Citation:**

Szalaj D, Silva A, Ré P and
Cabral H (2021) Detecting Regime
Shifts in the Portuguese Continental
Shelf Ecosystem Within the Last
Three Decades.
Front. Mar. Sci. 8:629130.
doi: 10.3389/fmars.2021.629130

<footnote n="1">¹ MARE – Marine and Environmental Sciences Centre, Faculdade de Ciências, Universidade de Lisboa, Lisbon, Portugal,</footnote>

<footnote n="2">² Instituto Português do Mar e da Atmosfera, Lisbon, Portugal, ³ Instituto Dom Luiz, Faculdade de Ciências, Universidade de Lisboa, Lisbon, Portugal, ⁴ INRAE, UR EABX, Centre Nouvelle-Aquitaine Bordeaux, Cestas, France</footnote>
# INTRODUCTION

Marine ecosystems worldwide are affected by increasing natural and anthropogenic pressures and consequently undergo significant changes at unprecedented rates. Affected by these changes, ecosystems can reorganise and still maintain the same function, structure, and identity (Walker et al., 2004). However, under some circumstances, the ecosystem may undergo changes that modify the system's structure and function and this process can be described as a shift to a new regime (Holling, 1973; May, 1977; Scheffer et al., 2001). Usually, a regime shift is triggered by large-scale climate-induced variations (Parsons and Lear, 2001), intense fishing exploitation (Daskalov, 2002) or both (Llope et al., 2011). Criteria used to define regime shifts vary and the changes that have to occur in order to consider that a system has undergone a regime shift are not well-defined (Lees et al., 2006). Normally, regime shifts are defined as high amplitude, low-frequency and often abrupt changes in species abundance and community composition that are observed at multiple trophic levels (TLs) (McKinnell et al., 2001). These changes are expected to occur on a large spatial scale and take place concurrently with physical changes in the climate system (McKinnell et al., 2001; Scheffer et al., 2001; Cury and Shannon, 2004; Collie et al., 2004; deYoung et al., 2004; Wooster and Zhang, 2004; Lees et al., 2006). Regime shifts have been described in several marine ecosystems including Northern Benguela (Heymans and Tomczak, 2016), the North Sea (Weijerman et al., 2005), and the Baltic Sea (Tomczak et al., 2013). In large upwelling ecosystems, it is common to observe decadal fluctuations in species abundance and their replacements (Cury et al., 2000). These fluctuations might be irreversible and might be an indicator of the new regime, as was the case in the Northern Benguela ecosystem (Heymans and Tomczak, 2016). However, changes in the upwelling systems might be interpreted as fluctuations within the limits of natural variability for an ecosystem, and not as an indicator of the regime shift (Cury and Shannon, 2004). The Portuguese continental shelf ecosystem (PCSE) constitutes the northernmost part of the Canary Current Upwelling System and is characterised by seasonal upwelling that occurs during the spring and summer as a result of steady northerly winds (Wooster et al., 1976; Fiúza et al., 1982). It has recently changed in the abundance of coastal pelagic species such as sardine, chub mackerel, horse mackerel, blue jack mackerel and anchovy (ICES, 2008, 2017a; Martins et al., 2013; Garrido et al., 2015). Moreover, in the last decades, an increase in higher TL species has been documented (Leitão, 2015a). The causes underlying changes in the pelagic community are not clear but it has been suggested that they result from a complex interplay between environmental variability, species interactions and fishing pressure (Leitão et al., 2014a; Leitão, 2015b; Veiga-Malta et al., 2019). There is evidence, that changes in the intensity of the Iberian coastal upwelling (resulting from the strengthening or weakening northern winds) had occurred in the last decades. However, the character of these changes is contradictory where some authors observed intensification of upwelling-favourable winds (Miranda et al., 2012; Pires et al., 2013) while others documented their weakening (Barton et al., 2013; Sydeman et al., 2014). In the latest review

of upwelling rate and intensity along the Portuguese coast, Leitão et al. (2019) documented a successive weakening of the upwelling since 1950, that lasted till mid/late 1970s in the north-west and south-west and till 1994 in the south coast. In recent years (1985–2009) an increase in upwelling index (UI) in all studied regions was documented while additionally upwelling intensification was observed in the south (Leitão et al., 2019). Moreover, a continuous increase in water temperature, ranging from 0.1 to 0.2°C decade⁻<fnref n="1" /> has also been documented (Baptista et al., 2017).

Understanding the changes that occurred in ecosystems, including their causes, timeline and impacts, is very important from the standpoint of ecosystem management as the potential of regime shifts raises several challenges for management (Crépin et al., 2012; Sguotti and Cormon, 2018). Therefore, identifying regime shifts and accounting for these changes are important in the context of Ecosystem-based management (Levin and Möllmann, 2015). This relates also to fisheries management because a regime shift might cause significant changes to the structure of the ecosystem and consequently impact the harvesting of fish stocks (Folke et al., 2004), as has been reflected in Portuguese fishery landings profiles in the last half-century (Leitão, 2015a).

Möllmann et al. (2009) used the integrated trend analysis (ITA) approach that combines multivariate statistics (principal component analysis, PCA) and methodologies to identify step changes in time-series such as sequential *t*-test analysis of regime shifts (STARS; Rodionov, 2004) and chronological clustering (CC; Legendre et al., 1985) to demonstrate the reorganization of the Baltic Sea ecosystem. By integrating hydroclimatic, nutrient, plankton and fisheries data into the ITA approach they identified two stable states between 1974 and 2005, separated by a transition regime in 1988–1993. Heymans and Tomczak (2016) and Tomczak et al. (2013) combined ITA with outputs of Ecopath with Ecosim modelling approach<fnref n="1" /> and used ecological network analysis (ENA) indicators as described by Ulanowicz (1986), to test for the changes in food web organisation and to describe changes that occurred in the Northern Benguela and the Baltic Sea.

The main objective of this study is to investigate the hypothesis that a regime shift occurred in the PCSE during the study period. Two types of data describing ecosystem drivers (fishing mortality and environment/climatic indices) and ecosystem state (biomass and ecosystem indicators), and ITA methodology (PCA, STARS, and CC) were used to: (1) assess the changes in the PCSE over the last three decades; (2) extract the common ecosystem trends (PCA); (3) detect significant shifts in datasets (STARS and CC); (4) compare times when the shifts occurred in different datasets; (5) explore the relationship between extracted trends, and (6) discuss the role of ecosystem drivers in the shifts describing ecosystem state.

<footnote n="1">¹www.ecopath.org</footnote>
# MATERIALS AND METHODS

## Study Site

The PCSE is situated in the northeast Atlantic between 36.5° and 42° N and between 10.5° and 7.5° W (Figure 1) and covers an area of around 25,000 km². Sea surface temperature (SST) in the Portuguese shelf waters varied from 13 and 17°C in the north and 15 and 21°C in the south in the winter and summer, respectively. The western Portuguese coast is affected by seasonal upwelling that occurs during spring and summer due to steady northerly winds (Wooster et al., 1976; Fiúza et al., 1982). Also, the Iberian Peninsula is influenced by mid-term atmospheric phenomena such as Atlantic multidecadal oscillation (AMO) and North Atlantic oscillation (NAO) (Trigo et al., 2004; Santos et al., 2011). In the Portuguese continental shelf, the most important fishing fleet in terms of volume of catch is the purse seine while the fleets that target more valuable species and are important due to the value of catch as opposed to volume, are bottom-trawl and multi-gear fleet. Purse-seine fisheries target mainly pelagic fish, where the most landed species are sardine (Sardina pilchardus), chub mackerel (Scomber colias), and horse mackerel (Trachurus trachurus) (Instituto Nacional de Estatística, 2000-2018). According to National Portuguese Statistical Institute (INE) statistics, bottom-trawlers that operate on the continental shelf target finfish while the multi-gear fleet targets sardine (but its catch is negligible), cephalopods, bivalves and crustaceans.

## Data and Model

In order to test the hypothesis of the occurrence of regime shifts in the PCSE during the last 30 years, we used data describing ecosystem drivers that include: (1) fishing mortality, (2) environmental variables and climatic indices, and data describing ecosystem state that includes: (1) observed biomass, (2) modelled biomass, and (3) ecosystem indicators (Table 1 and Figure 2).

Fishing mortality and observed biomass data used in this study were also used as inputs of the EwE model, describing the PCSE between 1986 and 2017. The static version of this model was developed by Veiga-Malta et al. (2019) and further updated and fitted to 32 years of time-series data by Szalaj et al. (2021). Moreover, modelled biomass and ecosystem indicators data used in the study come as an output of this EwE model. The model was developed using the EwE modelling software that uses a set of differential equations to describe the ecosystem's food web structure, its functioning, and dynamics. It has been described in detail by Christensen and Pauly (1992); Pauly et al. (2000), and Christensen and Walters (2004) (see text footnote 1). Using EwE data to assess the PCSE, allows for the integration of time-series data into the analysis that are not available from other sources, such as biomass of lower TL species and marine mammals. The EwE model estimated these groups and using them in the ITA provides a more complete picture of the studied ecosystem. In addition, the EwE model provided an extensive suite of indicators describing the ecosystem's properties, its structure and functioning. The use of these data allows for the assessment of the changes that occurred in the PCSE at the whole ecosystem level.

<figure><img src="image_3.png" /><figcaption>FIGURE 1 | Map showing limits of the study area (200 m isobath).</figcaption></figure>

Model description and details regarding model parametrization are available in the study by Szalaj et al. (2021). In summary, the model covers the area of the Portuguese Continental Shelf (0–200 m; ICES subdivisions IXa), consists of 33 functional groups (Figure 3) and includes three fishing fleets: bottom trawl (BT), pelagic and multi-gear. Functional groups aggregate the species that live in a similar environment and have a similar diet. A detailed description of species that comprise each functional group is available in the **Supplementary Table 1**.

## Ecosystem Drivers Data

Data describing ecosystem drivers consisted of two data groups: (1) fishing mortality and (2) environmental variables and climatic indices (Table 1). Fishing mortality (F) data for the majority of functional groups was calculated using Stock reduction analysis (SRA; Kimura et al., 1984) and implemented as a stochastic SRA (Walters et al., 2006) while F for sardine and hake was obtained from stock assessment outputs (ICES, 2018a,b). A detailed description of the SRA methodology and obtained results are provided in the **Supplementary Material S2**. Environmental variables consisted of annual, summer and winter means of SSTs and UI while climatic indices included annual, summer and winter means of AMO, NAO, and East Atlantic pattern (EA).
<table><caption>**TABLE 1** | Groups of data used in integrated trend analysis (ITA) and their sources.</caption><thead><tr><th>Data groups</th><th>Source</th></tr></thead><tbody><tr><td><strong>Ecosystem Drivers</strong></td><td></td></tr><tr><td><strong>Fishing mortality</strong></td><td></td></tr><tr><td>Demersal piscivorous fish, rays, benthopelagic piscivorous fish (BFP), sparids, bogue (Boops boops), horse mackerel (Trachurus trachurus), benthopelagic invertivorous fish (BFI), mackerel (Scomber scombrus), blue jack mackerel (Trachurus picturatus), and anchovy (Engraulis encrasicolus)</td><td>Stock reduction analysis (SRA) output</td></tr><tr><td>Sardine (Sardina pilchardus) and Hake (Merlucius merlucius)</td><td>ICES Stock assessment output (ICES, 2018a,b)<br>An estimate of fishing mortality (F)</td></tr><tr><td><strong>Environmental/climatic variables</strong></td><td></td></tr><tr><td><em>Environmental variables</em></td><td></td></tr><tr><td>Sea surface temperature (SST)</td><td>Met Office (https://www.metoffice.gov.uk/hadobs/hadsst3/data/download.html)</td></tr><tr><td>Upwelling index (UI)</td><td>Spanish Institute of Oceanography (IEO)<br>(http://www.indicedeafloramiento.ieo.es/interactivo.html)</td></tr><tr><td><em>Climatic indices</em></td><td></td></tr><tr><td>Atlantic multidecadal oscillation (AMO)</td><td>NOAA Earth System Research Laboratory<br>(http://www.esrl.noaa.gov/psd/data/timeseries/AMO)</td></tr><tr><td>North Atlantic oscillation (NAO)</td><td>NOAA Climate Prediction Centre-National Weather Service<br>(http://www.cpc.ncep.noaa.gov/data/teledoc/ea.shtml)</td></tr><tr><td>Eastern Atlantic pattern (EA)</td><td>NOAA Climate Prediction Centre-National Weather Service<br>(http://www.cpc.ncep.noaa.gov/data/teledoc/ea.shtml)</td></tr><tr><td><strong>ECOSYSTEM STATE</strong></td><td></td></tr><tr><td><strong>Observed biomass</strong></td><td></td></tr><tr><td><em>Demersal community</em></td><td></td></tr><tr><td>Hake (Merlucius merlucius),</td><td>Stock assessment output (ICES, 2018a)</td></tr><tr><td>Benthic cephalopods, sharks, demersal invertivorous fish, and flatfish</td><td>An estimate of ICES spawning stock biomass (SSB)</td></tr><tr><td>Rays, demersal piscivorous fish (DFP), Sparids</td><td>BT survey</td></tr><tr><td><em>Pelagic community</em></td><td>Stock reduction analysis output</td></tr><tr><td>Adult and juvenile sardine (Sardina pilchardus),</td><td>Stock assessment output (ICES, 2018b)</td></tr><tr><td>Horse mackerel (Trachurus trachurus)</td><td>For adult sardine and horse mackerel: an estimate of ICES spawning stock biomass (SSB)</td></tr><tr><td>Squids</td><td>For juvenile sardine estimate of ICES recruitment</td></tr><tr><td>Chub mackerel (Scomber colias)</td><td>BT survey</td></tr><tr><td>Bogue (Boops boops), mackerel (Scomber scombrus), blue jack mackerel (Trachurus picturatus), anchovy (Engraulis encrasicolus)</td><td>Acoustic survey</td></tr><tr><td><em>Benthopelagic community</em></td><td>Stock reduction analysis output</td></tr><tr><td>Benthopelagic piscivorous fish (BFP) and benthopelagic invertivorous fish (BFI)</td><td>Stock reduction analysis output</td></tr><tr><td><strong>Modelled biomass</strong></td><td></td></tr><tr><td><em>Mammals and birds</em></td><td></td></tr><tr><td>Bottlenose dolphin, Harbour porpoise, stripped dolphin, seabirds, common dolphin, and minke whale</td><td>EwE model output</td></tr><tr><td><em>Demersal community</em></td><td></td></tr><tr><td>Hake (Merlucius merlucius), rays, sharks, demersal invertivorous fish (DFI), demersal piscivorous fish (DFP), benthic cephalopods, sparids, and flatfish</td><td>EwE model output</td></tr><tr><td><em>Pelagic community</em></td><td></td></tr><tr><td>Tunas, squids, bogue (Boops boops), chub mackerel (Scomber colias), horse mackerel (Trachurus trachurus), mackerel (Scomber scombrus), blue jack mackerel (Trachurus picturatus), anchovy (Engraulis encrasicolus), juvenile sardine, and adult sardine (Sardina pilchardus)</td><td>EwE model output</td></tr><tr><td><em>Benthopelagic community</em></td><td></td></tr><tr><td>Benthopelagic piscivorous fish (BFP) and Benthopelagic invertivorous fish (BFI)</td><td>EwE model output</td></tr></tbody></table>

(Continued)
<table><caption>**TABLE 1** | Continued</caption><thead><tr><th>Data groups</th><th>Source</th></tr></thead><tbody><tr><td><strong>ECOSYSTEM STATE</strong></td><td></td></tr><tr><td><strong>Modelled biomass</strong></td><td></td></tr><tr><td>Invertebrates and Plankton</td><td></td></tr><tr><td>Henslow's crab (Polybius henslowi), shrimps, macrozoobenthos, suprabenthic invertebrates, macrozooplankton, meso/microzooplankton, and phytoplankton</td><td>EwE model output</td></tr><tr><td><strong>Ecosystem indicators</strong></td><td></td></tr><tr><td>Biomass-based</td><td></td></tr><tr><td>Demersal per pelagic biomass (Dem/Pel B), predatory biomass (Pred B), Kempton's index (Kempton Q), and Shannon diversity index (Shannon)</td><td>EwE model output</td></tr><tr><td>Catch-based</td><td></td></tr><tr><td>Total catch (Tot C), demersal per pelagic catch (Dem/Pel C), predatory catch (Pred C), pelagic catch (Pel C), and demersal catch (Dem C)</td><td>EwE model output</td></tr><tr><td>Trophic-based</td><td></td></tr><tr><td>Mean trophic level of the catch (mTLc), Mean trophic level of the community (mTLco)</td><td>EwE model output</td></tr><tr><td>Ecological network analysis (ENA)</td><td></td></tr><tr><td>Relative ascendency (A/C), Finn's cycling index (FCI), turnover rate (TotP/TotB), entropy (H), and redundancy (R)</td><td>EwE model output</td></tr></tbody></table>

For functional groups that consist of one species, a Latin name was provided in the brackets. For other functional groups, the list of species that are aggregated into the group is provided in **Supplementary Table 1**.

<figure><img src="image_4.png" /></figure>

## Ecosystem State Data

Data describing ecosystem state consisted of three data groups: (1) observed biomass, (2) modelled biomass, and (3) ecosystem indicators.

### Observed biomass data

Observed biomass data consisted of the biomass data aggregated into 20 functional groups (Table 1). These data were used to

calibrate the PCSE EwE model between 1986 and 2017 (Szalaj et al., 2021). When possible we used the outputs of stock assessment models as biomass estimates (i.e., sardine, hake and horse mackerel). For other species, when enough data were available, we estimated plausible biomass ($B_{SRA}$) trends using SRA. For the majority species (demersal piscivorous fish, rays, sparids, bogue, benthopelagic invertivorous, and piscivorous fish, mackerel, blue jack mackerel) we used BT survey data
<figure><img src="image_5.png" /></figure>

as abundance estimates for SRA input, with the exception of anchovy for which we used abundance estimates obtained from acoustic surveys. For species that lacked sufficient data to perform a SRA, the abundance index trends were taken directly from BT surveys (sharks, squids, benthic cephalopods, flatfish, and demersal invertivorous fish), with exception of the biomass for chub mackerel which was taken from an acoustic survey. Portuguese BT surveys are carried out in autumn (September–October) at predefined stations in the Portuguese Exclusive Economic Zone at depths varying from 20- to 500-m depth (ICES, 2017b). In this study only data that referred to the continental shelf (depth < 200 m) was considered. BT survey biomass estimates were calculated using the swept area method. Portuguese acoustic surveys (PELAGO survey series) are performed annually in spring (April–May) in an area that covers the Portuguese continental shelf. Acoustic fish data were obtained following the methodology described in ICES (2016) by performing echo-integration along 1 nm Elemental sampling distance unit (ESDU, distance measure used in acoustic surveys) and supported by fish sample collection. Both acoustic and BT surveys were performed by the Portuguese Institute of Sea and Atmosphere (IPMA).

### Modelled biomass data

Modelled biomass data for 33 functional groups that range from plankton to marine mammals (Figure 3 and Table 1) were obtained from the EwE model fitted to the time series for the PCSE between 1986 and 2017 (Szalaj et al., 2021). A full list of functional groups used in this study is presented in Supplementary Table 1.

### Ecosystem indicators data

Various indicators that describe ecosystem properties were extracted from the EwE model fitted to the time series for the PCSE between 1986 and 2017. Among the available indices, 16 indicators that represent distinctive processes, were used in the analysis. They are (1) demersal per pelagic biomass (Dem/Pel B), (2) predatory biomass (Pred B), (3) Kempton's index (Kempton Q), (4) Shannon diversity index (Shannon), (5) total catch (Tot C), (6) demersal per pelagic catch (Dem/Pel C), (7) predatory catch (Pred C), (8) pelagic catch (Pel C), (9) demersal catch (Dem C), (10) mean trophic level of the catch (mTLc), (11) mean trophic level of the community (mTLco), (12) relative ascendency (A/C), (13) Finn's cycling index (FCI), (14) turnover rate (TotP/TotB), (15) entropy (H), and (16) redundancy (R). Their definitions are presented in Table 2.

### Integrated Analysis

Data described in subsection “Data and Model” was subjected to an ITA that was performed using the method described by Diekmann and Möllmann (2010) (Figure 2). It integrates the following statistical methods: (1) PCA, (2) STARSs (Rodionov, 2004), and (3) CC (Legendre et al., 1985). The PCA and CC analyses were performed using R (R Core Team, 2020), packages vegan (Oksanen et al., 2019) and rioja (Juggins, 2020). The script was obtained from Diekmann et al. (2012). The STARS analysis was performed using STARS software written in visual basic for application (VBA) Excel 2002, available at: sites.google.com/climateologic/download.
<table><caption>TABLE 2 | Ecosystem indicators calculated from Ecopath with Ecosim (EwE), that describe the Portuguese continental shelf ecosystem (PCSE) between 1986 and 2017.</caption><thead><tr><th></th><th>Ecosystem indicators</th><th>Definition</th><th>Description</th></tr></thead><tbody><tr><td colspan="4"><strong>Biomass-based indicators</strong></td></tr><tr><td>1</td><td>Demersal per pelagic biomass (Dem/Pel B)</td><td>Sum of the biomass of demersal species/sum of the biomass of pelagic species</td><td>Explores processes benefiting the demersal or pelagic compartments of the ecosystem. In general, this ratio is expected to increase with fishing (Cury et al., 2005)</td></tr><tr><td>2</td><td>Predatory biomass (t-km<sup>-2</sup>·year<sup>-1</sup>) (Pred B)</td><td>Sum of the biomass of all species with a trophic level ≥ 4</td><td>Includes biomass of all the groups with TL ≥ 4 and tends to decrease with increasing fishing pressure in marine ecosystems (Rochet and Trenkel, 2003)</td></tr><tr><td>3</td><td>Kempton's index (Kempton Q)</td><td>Kempton Q = Q90 = <sup>0.8S</sup>&frasl;<sub>log(<sup>B<sub>2</sub></sup>&frasl;<sub>R<sub>1</sub></sub>)</sub></td><td>Where S is the total number of functional groups in the model; R1 and R2 are the representative biomass values of the 10th and 90th percentiles in the cumulative abundance distribution. It defines the inverse slope of the species-abundance curve<br>Describes diversity by considering functional groups with TL ≥ 3. This indicator tends to decrease with ecosystem degradation (Kempton and Taylor, 1976)</td></tr><tr><td>4</td><td>Shannon diversity index (Shannon)</td><td>Shannon = - &sum;<sub>i=1</sub><sup>N</sup> p<sub>i</sub>Lnp<sub>i</sub></td><td>Where N is the total number of species in the community ecosystem, p<sub>i</sub> is the proportion of total abundance represented by species i. The Shannon index increases as both the richness and the evenness of the community increase</td></tr><tr><td colspan="4"><strong>Catch-based indicators</strong></td></tr><tr><td>5</td><td>Total catch (t-km<sup>-2</sup>·year<sup>-1</sup>) (Tot C)</td><td>Tot C = &sum;<sub>i</sub> Y<sub>i</sub></td><td>Where Y<sub>i</sub> is the yield of compartment/species i<br>Includes total landings and discards of the different fleets and provides an idea of the overall fishing impact on the ecosystem (Pauly et al., 1998)</td></tr><tr><td>6</td><td>Demersal per pelagic catch (Dem/Pel C)</td><td>Sum of the catch of all demersal species/sum of the catch of all pelagic species</td><td>Explores fishing strategy and assess its focus on the demersal vs. pelagic compartment</td></tr><tr><td>7</td><td>Predatory catch (t-km<sup>-2</sup>·year<sup>-1</sup>) (Pred C)</td><td>Sum of the catch of all species with a trophic level ≥ 4</td><td>Expresses the catch of the groups with TL ≥ 4. It tends to increase with increasing fishing impact in the marine ecosystem</td></tr><tr><td>8</td><td>Pelagic catch (t-km<sup>-2</sup>·year<sup>-1</sup>) (Pel C)</td><td>Sum of the catch of all pelagic species</td><td>Includes total landings and discards of the fleets directed toward pelagic species. It provides an idea of the fishing impact on the pelagic ecosystem</td></tr><tr><td>9</td><td>Demersal catch (t-km<sup>-2</sup>·year<sup>-1</sup>) (Dem C)</td><td>Sum of the catch of all demersal species</td><td>Includes total landings and discards of the fleets directed toward demersal species. It provides an idea of the fishing impact on the demersal ecosystem</td></tr><tr><td colspan="4"><strong>Trophic-based indicators</strong></td></tr><tr><td>10</td><td>Mean trophic level of the catch (mTLc)</td><td>mTLc = <sup>&sum;<sub>i</sub> TL<sub>i</sub> &middot; Y<sub>i</sub></sup>&frasl;<sub>&sum;<sub>i</sub> Y<sub>i</sub></sub><br>TL<sub>j</sub> = 1 + &sum;<sub>j</sub> DC<sub>ij</sub> &middot; TL<sub>i</sub></td><td>Where Y<sub>i</sub> is the yield of compartment i; TL<sub>i</sub> is the trophic level of that compartment; TL<sub>j</sub> is the trophic level of predator j; TL<sub>i</sub> is a trophic level of prey i and DC<sub>ij</sub> is the proportion of prey i in the diet of predator j<br>Expresses the Trophic Level (TL) of the catch, reflects the fishing strategy of the fleet and is used to quantify the impact of fishing (Pauly et al., 1998). It is expected to decrease with the increase of fishing (Shannon et al., 2009)</td></tr><tr><td>11</td><td>Mean trophic level of the community (mTLco)</td><td>mTLco = &sum;<sub>i</sub> TL<sub>i</sub> <sup>B<sub>i</sub></sup>&frasl;<sub>B</sub></td><td>Where TL<sub>i</sub> is the trophic level of species B<sub>i</sub> is the species biomass and B is the total biomass<br>Describes the TL and reflects the structure of the ecosystem. It is used to quantify the impact of fishing (Rochet and Trenkel, 2003)</td></tr><tr><td colspan="4"><strong>Ecological network analysis (ENA) indicators</strong></td></tr><tr><td>12</td><td>Relative ascendency (%) (A/C)</td><td>A = &sum;<sub>i,j</sub> (T<sub>ij</sub>) &middot; log(<sup>T<sub>ij</sub> &middot; TST</sup>&frasl;<sub>T<sub>j</sub> &middot; T<sub>i</sub></sub>)<br>C = &sum;<sub>j</sub> T<sub>j</sub> log(<sup>T<sub>j</sub></sup>&frasl;<sub>TST</sub>)</td><td>Where A is ascendency where T<sub>ij</sub> is the flow between two compartments i and j and it includes all outflows from each compartment, T<sub>j</sub> is the sum of all material leaving the jth compartment, and T<sub>j</sub> is the sum of all flows entering the jth compartment. C is the development capacity where TST is the sum of all flows T<sub>ij</sub> in the system<br>Relative ascendency is an index if the organisation of the food web (Ulanowicz, 1986). It is negatively related to maturity (Christensen, 1995) which is related to the system's growth and development (Ulanowicz, 1986)</td></tr></tbody></table>

(Continued)
<table><caption>TABLE 2 | Continued</caption><thead><tr><th></th><th>Ecosystem indicators</th><th>Definition</th><th>Description</th></tr></thead><tbody><tr><td></td><td>Ecological network analysis (ENA) indicators</td><td></td><td></td></tr><tr><td>13</td><td>Finn's cycling index (%) (FCI)</td><td>FCI = <sup>TST<sub>c</sub></sup>&frasl;<sub>TST<sub>tot</sub></sub></td><td>Finn's cycling index is the proportion of the total system throughput (TST<sub>tot</sub>) that is recycled in the system/TST<sub>c</sub>). Quantifies the relative amount of recycling in the system and is an indicator of stress and structural differences (Finn, 1976). It is considered an indicator of ecosystem's ability to maintain its structure and integrity (Monaco and Ulanowicz, 1997) and an indicator of stress (Ulanowicz, 1986). An increase in FCI would suggest that the system has a capability to recover faster from a perturbation, due to an increased system's recycling and consequently higher capacity to conserve nutrients. It is important to keep in mind that a system would be expected to take longer to recover (lower FCI) when it is in a more degraded state (Vasconcellos et al., 1997)</td></tr><tr><td>14</td><td>Turnover rate (TotP/TotB)</td><td>The total production/total biomass ratio</td><td>System turnover rate is an indicator of the average size of organisms in an ecosystem (Shannon et al., 2009). It is hypothesized that fishing reduces the mean size of organisms in ecosystems therefore the turnover rate is expected to increase with fishing (Shannon et al., 2009)</td></tr><tr><td>15</td><td>Entropy (H)</td><td>H = - &sum;<sub>j</sub> <sup>T<sub>ij</sub></sup>&frasl;<sub>TST</sub> &middot; log(<sup>T<sub>ij</sub></sup>&frasl;<sub>TST</sub>)</td><td>Where T<sub>ij</sub> is the flow between two compartments and TST is the sum of all flows in the system<br>Entropy represents the total number and diversity of flows in the system (Mageau et al., 1998). The diversity of flows is similar to diversity assessed by the Shannon diversity index (Shannon-Wiener) but considered flows instead of abundance (Ulanowicz and Norden, 1990)</td></tr><tr><td>16</td><td>Redundancy (R)</td><td>R = - &sum;<sub>i=1</sub><sup>n</sup> &sum;<sub>j=1</sub><sup>n</sup> (T<sub>ij</sub> &middot; log(<sup>T<sub>ij</sub><sup>2</sup></sup>&frasl;<sub>&sum;<sub>i=1</sub><sup>n</sup> T<sub>ij</sub> &middot; &sum;<sub>i=1</sub><sup>n</sup> T<sub>ij</sub></sub>)</td><td>Where: T<sub>ij</sub> is the flow between two compartments<br>This index is considered as an index of the system's resilience (Heymans et al., 2007) and it has been linked with ecosystem stability (Christensen, 1995). It describes the distribution of energy flow pathways in the system. If the R is high, it indicates that they are many alternative pathways for energy flows to get from one compartment/functional group to another that allows compensation for the environmental stress, while lower R indicate that flows are constrained and energy is channelled through more specific pathways</td></tr></tbody></table>

The PCA, based on the correlation matrix, was performed separately on fishing mortality data (PCA_FM), environmental/climatic variables (PCA_ENV_CLIM), observed biomass (PCA_B_objs), modelled biomass (PCA_B_model) and ecosystem indicators (PCA_indexes). Before analysis, each dataset was checked for cross-correlation (Pearson) and variables that are highly correlated (>0.8, Supplementary Material S3) were excluded from the PCA. The list of the variables included in the PCA analysis is shown in Table 3. The number of meaningful principal components (PC) that were considered further in the analysis was assessed by the *scree plot* that shows how much variation each PC captures from the data and the *broken stick plot* that demonstrates how many PCs should be retained. Also, variable loadings for each meaningful PC were checked and the variables whose loadings were < 0.3 were excluded from further analysis. Year scores for the meaningful principal components (mainly PC1 and PC2) were plotted against time to display temporal relationships (PC trends). Then PC trends were subjected to STARS analysis to detect the occurrence of regime shifts in the ecosystem. Moreover, the relationship between PC trends extracted from the data describing ecosystem state and ecosystem drivers were explored with cross-correlation. To identify significant shifts in mean values of a given time

series, STARS was performed separately on each time series of the whole dataset and on the time series of PC scores (PC trends). The STARS method is based on a sequential analysis of change points in the time series using the Student t-test. At each observation, the procedure calculates a regime shift index (RSI) to determine whether the following observation values are significantly different from the mean of the previous regime. The method and its use in the detection of shifts in marine systems were discussed in Strinimann et al. (2019). In order to account for serial correlations (red noise component), the prewhitening procedure that adjusts the significance level of the shifts by calculating effective degrees of freedom, described in Rodionov (2006), was used. In this study the red noise component was modelled by the first-order autoregressive model (AR1) and estimated based on the ordinary least squares (OLS) as suggested by Marriot and Pope (1954) and Kendall (1954). Parameters that control the scale and magnitude of potential regime shifts (significance level and cut-off length) and weights assigned to outliers (Huber's weight parameter) were set a priori to 0.05, 10, and 3, respectively according to Tomczak et al. (2013) and Heymans and Tomczak (2016). The magnitude of the detected regime shifts was assessed using RSI when the positive value indicates an increase and the negative value indicates a decrease.
<table><caption>TABLE 3 | Variables used in principal component analysis (PCA) after selection process (criteria for selection: cross-correlation ≤ 0.8).</caption><thead><tr><th rowspan="3">Name of PCA analysis<br>Type of data used in PCA analysis</th><th colspan="2">Ecosystem drivers</th><th colspan="3">Ecosystem state</th></tr><tr><th>Fishing mortality</th><th>Environmental/climatic variables</th><th>Observed biomass</th><th>Modelled biomass</th><th>Ecosystem indicators</th></tr></thead><tbody><tr><td rowspan="14">Variables used in PCA analysis (Cor ≤ 0.8)</td><td>Sardine</td><td>SSTm</td><td>Adult sardine</td><td>Adult sardine</td><td>Dem/Pel B</td></tr><tr><td>Anchovy</td><td>SSTw</td><td>Juvenile sardine</td><td>Juvenile sardine</td><td>Shannon</td></tr><tr><td>Horse mackerel</td><td>AMOm</td><td>Anchovy</td><td>Anchovy</td><td>Tot C</td></tr><tr><td>Mackerel</td><td>NAOm</td><td>Chub mackerel</td><td>Horse mackerel</td><td>Pred C</td></tr><tr><td>Blue jack mackerel</td><td>NAOw</td><td>Mackerel</td><td>Hake</td><td>TL catch</td></tr><tr><td>Hake</td><td>NAOs</td><td>Blue jack mackerel</td><td>Rays</td><td>R</td></tr><tr><td>Rays</td><td>EAm</td><td>Tunas</td><td>Benthopelagic invertivorous fish</td><td>FCI</td></tr><tr><td>Benthopelagic piscivorous fish</td><td>EAw</td><td>Hake</td><td>Benthic Cephalopods</td><td>H</td></tr><tr><td>Benthopelagic invertivorous fish</td><td>EAs</td><td>Rays</td><td>Squids</td><td>TotP/TotB</td></tr><tr><td rowspan="6">Demersal piscivorous fish</td><td>Ulm</td><td></td><td>Benthopelagic piscivorous fish</td><td>Stripped dolphin</td><td></td></tr><tr><td>Ulw</td><td></td><td>Benthopelagic invertivorous fish</td><td>Suprabenthic invertrebrates</td><td></td></tr><tr><td>Uls</td><td></td><td>Demersal invertivorous fish</td><td></td><td></td></tr><tr><td></td><td></td><td>Sparids</td><td></td><td></td></tr><tr><td></td><td></td><td>Flatfish</td><td></td><td></td></tr><tr><td></td><td></td><td>Sharks</td><td></td><td></td></tr><tr><td></td><td></td><td>Benthic cephalopods</td><td></td><td></td></tr><tr><td></td><td></td><td>Squids</td><td></td><td></td></tr><tr><td>Number of variables used</td><td>10</td><td>12</td><td>17</td><td>11</td><td>9</td></tr></tbody></table>

The CC discontinuity analysis is a clustering technique that groups only sequential samples and has been used to identify the abrupt changes in time-series. It is used for multivariate datasets and clustering is based on the calculation of a similarity matrix of the time-series based on Euclidian distances. The CC analysis was performed separately on the same data groups as the PCA (fishing mortality, environmental/climatic variables, observed biomass, modelled biomass, and ecosystem indicators). According to Diekmann and Möllmann (2010), the significance level α (considered as a clustering-intensity parameter) and connectedness level were set to 0.01 and 50%, respectively. To make results interpretable in relation to PCA, data were standardised prior to the calculation of the similarity matrix.

# RESULTS

## Ecosystem Changes, Key Functional Groups, and Key Drivers (PCA)

In PCA_FM that included fishing mortality variables, three PCs that in total explained 72.8% of variance were retained (Table 4). The time trajectory of PC1_FM presented a relatively stable

trend until around 1998, followed by a decrease till the end of the study period (Figure 4A). Variables that contributed the most and positively related to this PC trend were fishing mortality of demersal piscivorous fish and rays (Table 4). The trend of PC2_FM showed an initial decrease until 1996, followed by an increase till 2006 and a further slight decrease until the end of the study period (Figure 4A). The variables that positively contributed to this trend were fishing mortality of benthopelagic invertivorous fish and mackerel, while fishing mortality of benthopelagic piscivorous fish, anchovy and hake contributed negatively (Table 4).

The trend of PC3_FM showed an initial gradual decrease with time till about 2008, followed by a sharp increase till the end of the study period (Figure 4B). The variables that contributed positively to this PC trend were fishing mortality of blue jack mackerel, horse mackerel (both positively) and sardine (negatively) (Table 4).

In the PCA_ENV_CLIM that included environmental/climatic variables, two retained PCs explained in total 54.6% of the variance (Table 4). The trend of PC1_ENV_CLIM showed a gradual increase with time and variables SST, AMO and EA showed the largest positive
<table><caption>TABLE 4 | Contribution of each variable to meaningful principal components for each PCA performed.</caption><thead><tr><th>Variable abbreviation</th><th>Variable name</th><th colspan="4">Meaningful PCs loadings</th></tr><tr><th>PCA_FM</th><th>PCA on fishing mortality</th><th>PC1</th><th>PC2</th><th>PC3</th><th>Total</th></tr></thead><tbody><tr><td>Variance explained (%)</td><td></td><td><strong>33.8</strong></td><td><strong>22.7</strong></td><td><strong>16.2</strong></td><td><strong>72.8</strong></td></tr><tr><td>F_DFP</td><td>Demersal piscivorous fish</td><td><strong>1.2</strong></td><td>-0.4</td><td>0.1</td><td></td></tr><tr><td>F_Rays</td><td>Rays</td><td><strong>1.3</strong></td><td>0.2</td><td>0.0</td><td></td></tr><tr><td>F_BFP</td><td>Benthopelagic piscivorous fish</td><td><strong>1.0</strong></td><td>-0.8</td><td>-0.2</td><td></td></tr><tr><td>F_HOM</td><td>Horse mackerel</td><td>0.8</td><td>0.4</td><td><strong>0.8</strong></td><td></td></tr><tr><td>F_BFI</td><td>Benthopelagic invertivorous fish</td><td>0.2</td><td><strong>0.8</strong></td><td>-0.6</td><td></td></tr><tr><td>F_MAC</td><td>mackerel</td><td>0.8</td><td><strong>0.9</strong></td><td>0.3</td><td></td></tr><tr><td>F_JAA</td><td>Blue jack mackerel</td><td>-0.6</td><td>-0.6</td><td><strong>0.7</strong></td><td></td></tr><tr><td>F_ANE</td><td>Anchovy</td><td>0.2</td><td><strong>-0.8</strong></td><td>0.7</td><td></td></tr><tr><td>F_PIL</td><td>Sardine</td><td>0.3</td><td>0.0</td><td><strong>-0.5</strong></td><td></td></tr><tr><td>F_HKE</td><td>Hake</td><td>0.4</td><td><strong>-0.8</strong></td><td>-0.8</td><td></td></tr><tr><td>PCA_ENV_CLIM</td><td>PCA on environmental/climatic variables</td><th>PC1</th><th>PC2</th><th></th><th></th></tr><tr><td>Variance explained (%)</td><td></td><td><strong>31.2</strong></td><td><strong>23.4</strong></td><td></td><td><strong>54.6</strong></td></tr><tr><td>AMOm</td><td>Annual average of Atlantic multidecadal oscillation</td><td><strong>1.0</strong></td><td>-0.2</td><td></td><td></td></tr><tr><td>NAOm</td><td>Annual average of North Atlantic oscillation</td><td>-0.7</td><td><strong>0.9</strong></td><td></td><td></td></tr><tr><td>NAOw</td><td>Winter average of North Atlantic oscillation</td><td>-0.3</td><td><strong>1.1</strong></td><td></td><td></td></tr><tr><td>NAOs</td><td>Summer average of North Atlantic oscillation</td><td><strong>-0.5</strong></td><td>-0.1</td><td></td><td></td></tr><tr><td>EAm</td><td>Annual average of Eastern Atlantic pattern</td><td><strong>1.0</strong></td><td>0.3</td><td></td><td></td></tr><tr><td>EAw</td><td>Winter average of Eastern Atlantic pattern</td><td><strong>0.5</strong></td><td>0.2</td><td></td><td></td></tr><tr><td>EAs</td><td></td><td><strong>1.0</strong></td><td>0.3</td><td></td><td></td></tr><tr><td>SSTm</td><td>Annual average sea surface temperature</td><td><strong>0.9</strong></td><td>0.7</td><td></td><td></td></tr><tr><td>SSTw</td><td>Winter average sea surface temperature</td><td>0.6</td><td><strong>0.8</strong></td><td></td><td></td></tr><tr><td>UI_m</td><td>Annual average upwelling index</td><td>-0.2</td><td><strong>0.7</strong></td><td></td><td></td></tr><tr><td>UI_w</td><td>Winter average upwelling index</td><td>-0.6</td><td><strong>0.9</strong></td><td></td><td></td></tr><tr><td>UI_s</td><td>Summer average upwelling index</td><td><strong>0.7</strong></td><td>-0.2</td><td></td><td></td></tr><tr><td>PCA_B_objs</td><td>PCA on observed biomass</td><th>PC1</th><th>PC2</th><th>PC3</th><th></th></tr><tr><td>Variance explained (%)</td><td></td><td><strong>34.1</strong></td><td><strong>23.1</strong></td><td><strong>16.1</strong></td><td><strong>80.7</strong></td></tr><tr><td>B_ANE_o</td><td>Anchovy</td><td><strong>-0.8</strong></td><td>0.6</td><td>-0.1</td><td></td></tr><tr><td>B_PIL_o</td><td>Adult sardine</td><td><strong>1.0</strong></td><td>-0.1</td><td>-0.1</td><td></td></tr><tr><td>B_PIL_Juv_o</td><td>Juvenile sardine</td><td><strong>1.0</strong></td><td>0.0</td><td>0.0</td><td></td></tr><tr><td>B_MAC_o</td><td>Mackerel</td><td><strong>-1.1</strong></td><td>-0.2</td><td>0.0</td><td></td></tr><tr><td>B_MASo</td><td>Chub mackerel</td><td><strong>-0.5</strong></td><td>-0.5</td><td>-0.3</td><td></td></tr><tr><td>B_Sparids_o</td><td>Sparids</td><td><strong>-1.1</strong></td><td>-0.2</td><td>0.0</td><td></td></tr><tr><td>B_Flatfish_o</td><td>Flatfish</td><td>0.1</td><td>0.2</td><td><strong>0.9</strong></td><td></td></tr><tr><td>B_Hake_o</td><td>Hake</td><td>0.1</td><td><strong>1.1</strong></td><td>-0.2</td><td></td></tr><tr><td>B_JAA_o</td><td>Blue jack mackerel</td><td>-0.5</td><td><strong>0.7</strong></td><td>-0.3</td><td></td></tr><tr><td>B_DFI_o</td><td>Demersal invertivorous fish</td><td>-0.2</td><td>0.4</td><td><strong>1.0</strong></td><td></td></tr><tr><td>B_BFI_o</td><td>Benhopelagic invertivorous fish</td><td><strong>0.8</strong></td><td>0.3</td><td>0.3</td><td></td></tr><tr><td>B_BFP_o</td><td>Benthopelagic piscivorous fish</td><td>-0.2</td><td><strong>1.0</strong></td><td>-0.1</td><td></td></tr><tr><td>B_Ben_Ceph_o</td><td>Benthic cephalopods</td><td>0.0</td><td>-0.2</td><td><strong>0.9</strong></td><td></td></tr><tr><td>B_Squids_o</td><td>Squids</td><td>0.4</td><td><strong>0.9</strong></td><td>0.1</td><td></td></tr><tr><td>B_Rays_o</td><td>Rays</td><td><strong>-1.0</strong></td><td>0.5</td><td>0.1</td><td></td></tr><tr><td>B_Sharks_o</td><td>Sharks</td><td>-0.2</td><td>0.1</td><td><strong>0.9</strong></td><td></td></tr><tr><td>B_Tunas_o</td><td>Tunas</td><td><strong>0.8</strong></td><td>0.7</td><td>-0.2</td><td></td></tr><tr><td>PCA_B_model</td><td>PCA on modelled biomass</td><th>PC1</th><th>PC2</th><th></th><th></th></tr><tr><td>Variance explained (%)</td><td></td><td><strong>47.0</strong></td><td><strong>32.4</strong></td><td></td><td><strong>79.4</strong></td></tr><tr><td>B_Supra_inv_m</td><td>Suprabenthic invertebrates</td><td>-0.3</td><td><strong>-0.9</strong></td><td></td><td></td></tr><tr><td>B_PIL_m</td><td>Adult sardine</td><td><strong>0.9</strong></td><td>-0.7</td><td></td><td></td></tr></tbody></table>

(Continued)
<table><caption>TABLE 4 | Continued</caption><thead><tr><th>Variable abbreviation</th><th>Variable name</th><th colspan="2">Meaningful PCs loadings</th></tr></thead><tbody><tr><td>PCA_B_model</td><td>PCA on modelled biomass</td><td>PC1</td><td>PC2</td></tr><tr><td>B_PIL_Juv_m</td><td>Juvenile sardine</td><td><strong>0.9</strong></td><td>-0.7</td></tr><tr><td>B_ANE_m</td><td>Anchovy</td><td><strong>-1.0</strong></td><td>-0.4</td></tr><tr><td>B_BFI_m</td><td>Benthopelagic invertivorous fish</td><td>0.1</td><td><strong>-1.1</strong></td></tr><tr><td>B_HOM_m</td><td>Horse mackerel</td><td><strong>-1.3</strong></td><td>0.0</td></tr><tr><td>B_Ben_Ceph_m</td><td>Benthic cephalopods</td><td><strong>-1.1</strong></td><td>0.1</td></tr><tr><td>B_Squids_m</td><td>Squids</td><td><strong>-1.2</strong></td><td>0.1</td></tr><tr><td>B_Rays_m</td><td>Rays</td><td><strong>-1.1</strong></td><td>-0.6</td></tr><tr><td>B_Hake_m</td><td>Hake</td><td>0.1</td><td><strong>-1.2</strong></td></tr><tr><td>B_Stripped_dolphin_m</td><td>Stripped dolphin</td><td>-0.5</td><td><strong>-0.9</strong></td></tr><tr><td>PCA_indices</td><td>PCA on ecosystem indicators</td><td>PC1</td><td>PC2</td></tr><tr><td>Variance explained (%)</td><td></td><td><strong>51.0</strong></td><td><strong>33.7</strong></td></tr><tr><td>Dem_Pel_B</td><td>Demersal per Pelagic biomass</td><td><strong>-1.3</strong></td><td>0.2</td></tr><tr><td>Shannon</td><td>Shannon diversity index</td><td><strong>-1.0</strong></td><td>0.1</td></tr><tr><td>Tot C</td><td>Total catch</td><td><strong>-1.0</strong></td><td>0.8</td></tr><tr><td>Pred C</td><td>Predatory catch</td><td>-0.6</td><td><strong>1.1</strong></td></tr><tr><td>mTLc</td><td>The trophic level of the catch</td><td><strong>1.1</strong></td><td>0.0</td></tr><tr><td>R</td><td>Redundancy</td><td>-0.7</td><td><strong>-1.2</strong></td></tr><tr><td>FCI</td><td>Fins cycling index</td><td><strong>-1.2</strong></td><td>-0.5</td></tr><tr><td>H</td><td>Entropy</td><td>-0.2</td><td><strong>1.3</strong></td></tr><tr><td>TotP_TotB</td><td>Turnover rate</td><td><strong>-1.1</strong></td><td>-0.6</td></tr></tbody></table>

Loadings values in bold showed the highest contribution.

contribution to this trend (Table 4 and Figure 4C). The time trajectory of PC2_ENV_CLIM fluctuated till 2010 without any clear trend, then increased substantially (Figure 4C). This trend was mainly explained by NAO and UI (both annual and winter averages) (Table 4).

In the PCA_Bobs that included EwE input biomass data, three retained PCs explained 73.5% of the variance (Table 4). The time trajectory of PC1_Bobs which explained 33.8% of the variance, showed a decreasing trend that represented the biomass trend of sardine and benthopelagic invertivorous fish (Table 4 and Figure 4D). The variables that showed an opposite trend to PC1_Bobs were pelagic fish (anchovy and mackerels), sparids and rays (Table 4). Time trajectory of PC2_Bobs that explained 23.7% of the variance and showed an initial sharp decrease until 1998 followed by a decrease (Figure 4D) represented the biomass trend for blue jack mackerel, hake, demersal piscivorous fish and squids (Table 4). The trend of PC3_Bobs that explained 15.9% of the variance, did not show any clear tendency and reflected the biomass of demersal invertivorous fish and benthic cephalopods (Table 4 and Figure 4E).

In the PCA_B_model that included EwE output biomass data, two retained PCs explained 79.4% of the variance. The time trajectory of PC1_B_model fluctuated without any clear trend till around 2006, when it substantially decreased (Figure 4F). It represented the modelled biomass trend for sardine (Table 4). The opposite trend was observed for other pelagics (horse mackerel and anchovy), cephalopods and rays. The trend of PC2_B_model increased till 1996, then stabilised and fluctuated till the end of the study period (Figure 4F).

It reflected the opposite trend to modelled biomass of hake, benthopelagic invertivorous fish, suprabenthic invertebrates and striped dolphin (Table 4).

Both PC1_indices and PC2_indices explained 51.0 and 33.7% of the total variance, respectively (Table 4). The time trajectory of PC1_indices showed an increasing trend (Figure 4G) that reflected the trend of mTLc (Table 4). The opposite trend to PC1_indices was observed for Tot C, Dem/Pel B, diversity (Shannon), and ENA indices (FCI, and TotP/TotB) (Table 4). The trend of PC2_indices sharply decreased until 1996, stabilised and since 2010 increased again (Figure 4G). The variables that reflected this trend were Pred C and H while R showed the opposite trend (Table 4).

## Relationship Between Ecosystem Drivers and Ecosystem State

The PC trends extracted from the observed biomass data group showed significant correlations with PC trends extracted from fishing mortality and environmental/climatic data groups (Table 5). The statistical correlation was observed between PC1_Bobs and PC1_FM (R = 0.82, *p* < 0.01), PC1_ENV_CLIM (R = -0.66, *p* < 0.01) and PC3_FM (R = 0.35, *p* = 0.05). The largest correlations for PC2_Bobs and PC3_Bobs were found for PC2_FM (R = 0.68, *p* < 0.01) and PC2_ENV_CLIM (R = 0.38, *p* = 0.03), respectively.

A statistical correlation was also found between PC trends extracted from modelled biomass and ecosystem drivers' data (Table 5). Here, the PC1_B_model had the largest statistically
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

**FIGURE 4** | Results of principal component analysis (PCA); time-trajectory of PC scores of PCA performed on (A) fishing mortality PC1 and PC2 (PCA_FM_PC1_PC2); (B) fishing mortality PC3 (PCA_FM_PC3); (C) environmental/climatic variables PC1 and PC2 (PCA_ENV_CLIM); (D) observed biomass PC1 and PC2 (PCA_B_obser PC1_PC2); (E) observed biomass PC3 (PCA_B_obser PC3); (F) modelled biomass PC1 and PC2 (PCA_B_model); and (G) ecosystem indicators PC1 and PC2 (PCA_indices) (black dots – PC1, PC3 white dots – PC2).
<table><caption>**TABLE 5** | Spearman cross-correlation between PCs extracted from different groups of data.</caption><thead><tr><th></th><th>PC1_B</th><th>PC2_B</th><th>PC3_B</th><th>PC1_B</th><th>PC2_B</th><th>PC1</th><th>PC2</th><th>PC3</th><th>PC1_ENV</th><th>PC2_ENV</th><th>PC1</th><th>PC2</th></tr><tr><th></th><th>_obs</th><th>_obs</th><th>_obs</th><th>_model</th><th>_model</th><th>_FM</th><th>_FM</th><th>_FM</th><th>_CLIM</th><th>_CLIM</th><th>_indices</th><th>_indices</th></tr></thead><tbody><tr><td>PC1_B obs</td><td>1</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PC2_B obs</td><td>0.02</td><td>1</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PC3_B obs</td><td>-0.14</td><td>0.10</td><td>1</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PC1_B model</td><td><strong>0.55</strong></td><td><strong>-0.77</strong></td><td>-0.06</td><td>1</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PC2_B model</td><td><strong>-0.44</strong></td><td><strong>-0.61</strong></td><td>-0.20</td><td>0.19</td><td>1</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PC1_FM</td><td><strong>0.82</strong></td><td>-0.33</td><td>-0.20</td><td><strong>0.76</strong></td><td>-0.2</td><td>1</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PC2_FM</td><td>0.32</td><td><strong>0.68</strong></td><td>-0.09</td><td><strong>-0.38</strong></td><td><strong>-0.37</strong></td><td>-0.04</td><td>1</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PC3_FM</td><td><strong>0.35</strong></td><td>0.33</td><td>0.22</td><td>-0.02</td><td><strong>-0.78</strong></td><td>0.21</td><td>0.01</td><td>1</td><td></td><td></td><td></td><td></td></tr><tr><td>PC1_ENV_CLIM</td><td><strong>-0.66</strong></td><td>-0.2</td><td>-0.03</td><td>-0.32</td><td><strong>0.43</strong></td><td><strong>-0.62</strong></td><td><strong>-0.40</strong></td><td>-0.16</td><td>1</td><td></td><td></td><td></td></tr><tr><td>PC2_ENV_CLIM</td><td>-0.29</td><td>0.31</td><td><strong>0.38</strong></td><td><strong>-0.46</strong></td><td>-0.27</td><td><strong>-0.41</strong></td><td>0.02</td><td>0.27</td><td>0.02</td><td>1</td><td></td><td></td></tr><tr><td>PC1_indices</td><td><strong>-0.88</strong></td><td>0.04</td><td>0.26</td><td><strong>-0.52</strong></td><td><strong>0.34</strong></td><td><strong>-0.86</strong></td><td>-0.22</td><td>-0.25</td><td><strong>0.62</strong></td><td>0.27</td><td>1</td><td></td></tr><tr><td>PC2_indices</td><td>0.14</td><td><strong>0.86</strong></td><td>0.03</td><td><strong>-0.62</strong></td><td><strong>-0.70</strong></td><td>-0.11</td><td><strong>0.55</strong></td><td><strong>0.54</strong></td><td>-0.23</td><td><strong>0.41</strong></td><td>-0.14</td><td>1</td></tr></tbody></table>

Values in bold indicate significant correlations ($p \le 0.05$).

significant correlation with PC1_FM ($R = 0.76, p < 0.01$) followed by PC2_ENV_CLIM ($R = -0.46; p < 0.01$) and PC2_FM ($R = -0.38, p = 0.03$). On the other hand, PC2_B_model shown the largest correlation with PC3_FM ($R = -0.78, p < 0.01$), followed by PC1_ENV_CLIM ($R = 0.43, p = 0.01$) and PC2_FM ($R = -0.37, p = 0.04$).

Also, statistical correlations were found between the PC trends extracted from the ecosystem indicators data group and both types of ecosystem drivers (Table 5). PC1_indices had the largest correlation with PC1_FM ($R = -0.86, p < 0.01$) followed by PC1_ENV_CLIM ($R = 0.62, p < 0.01$) while PC2_indices showed the largest correlation with PC2_FM ($R = 0.55, p < 0.01$), PC3_FM ($R = 0.54, p < 0.01$) followed by PC2_ENV_CLIM ($R = 0.41, p = 0.02$).

## Regime Shift Detection (STARS and CC)

Significant shifts detected by STARS were observed for many variables representing ecosystem drivers and ecosystem state (Figures 5–8 and Supplementary Table 4).

Among fishing mortality variables, significant shifts in the mid-90s were observed for fishing mortality of sparids, bogue, mackerel, benthopelagic piscivorous fish, and anchovy (Figure 5). The second shift in the mid-2000s was observed for fishing mortality of piscivorous fish, rays, sparids, and among small pelagic fish (SPF) for horse mackerel and anchovy (Figure 5). The last shift in the mid-2010s was observed for fishing mortality of all species (excluding fishing mortality of piscivorous fish and bogue).

Considering environmental/climatic variables, STARS detected significant shifts in the SST, AMO, EA in the mid-90s, followed by mid-2010s (2011–2016) (Figure 5). For NAO and UI, STARS detected only one shift in the mid-2010s (2014–2016) (Figure 5).

Among observed biomass variables, STARS detected significant shifts in the mid-90s for sardine, anchovy, mackerel, sparids, piscivorous fish, squids and hake. The second shift in the mid-2000s was observed for benthopelagic invertivorous fish, bogue, sparids, and rays. The last shift in the mid-2010s was

observed for benthopelagic invertivorous fish and SPF (horse mackerel and bogue) (Figure 6). For modelled biomass variables, STARS detected significant shifts mainly in the mid-2000s and mid-2010s. In the mid-90s a shift in biomass was detected only for sardine, anchovy, blue jack mackerel, bogue, and benthopelagic piscivorous fish. In the mid-2000s the shifts were detected in the biomass of cephalopods, lower TL groups and invertebrates, chub mackerel and fish feeding on invertebrates. Further shifts in the mid-2010s, were detected for sardine and other SPF (horse mackerel, chub mackerel, and bogue), cephalopods and lower TL functional groups (Figure 7).

For ecosystem indices, significant shifts were mainly detected in the mid-90s and mid-2010s. The shift in the mid-2000s was only detected for Dem/Pel C and Kempton Q (Figure 8). Shifts in the mid-90s were detected for catch related indices (Tot C, Dem C, Pred C, and Dem/Pel C), Dem/Pel B, Kempton Q and H (Figure 8). Shifts after 2010 were detected for Pred B, diversity indices (Kempton Q, Shannon), catch related indices (Dem C, Pel C, and mTLc) and ENA indices (FCI, R, H, and Tot P/Tot B) (Figure 8).

Shifts at similar times were also detected by STARS performed on PCA results (PC trends) and CC analysis performed on separate data groups (Table 6). When performed on PC trends extracted from the fishing mortality data, STARS indicated that significant shifts occurred mainly in 1995 and 1999 (for PC2_FM and PC1_FM, respectively), 2003 (for PC2_FM), and 2012–2016 (for PC1_FM and PC2_FM), while CC detected discontinuity in the fishing mortality data in 1992 and 2003 (Table 5).

For STARS performed on PC trends of environmental/climatic variables significant shifts were detected in 1997 for PC1_ENV_CLIM and 2014–2016 for PC2_ENV_CLIM and PC1_ENV_CLIM, respectively while CC detected discontinuity in 1995, 2002, and 2009 (Table 6).

Considering PC trends of observed biomass, STARS detected significant shifts in 1993 for PC2_B_objs, 2003–2005 for PC2_B_objs and PC1_B_objs, respectively and 2013–2015 for PC1_B_objs and PC2_B_objs (Table 6). The CC detected significant discontinuity only in 1990 (Table 6).
<figure><img src="image_13.png" /><figcaption>**FIGURE 5** | Time-series of ecosystem drivers data that includes fishing mortality and environmental/climatic variables. Significant shifts detected by STARS analysis are displayed on the plots as solid and dashed vertical line for shifts significant at *p* ≤ 0.05 and ≤0.1, respectively. Abbreviations, used in the figure, are explained in **Supplementary Table 4**.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**FIGURE 6** | Time-series observed biomass data. Significant shifts detected by STARS analysis are displayed on the plots as solid and dashed vertical line for shifts significant at *p* ≤ 0.05 and ≤0.1, respectively. Abbreviations, used in the figure, are explained in **Supplementary Table 4**.</figcaption></figure>

When performed on PC trends extracted from modelled biomass data group, STARS detected significant shifts in 1992 (for PC2_B_model), 2008 (for PC1_B_model), and 2015 (for PC1_B_model) while, CC showed discontinuity in the modelled biomass data in 1998 and 2013 (Table 6).

For STARS performed on PC trends of ecosystem indicators, significant shifts in the ecosystem were detected in 1992 for PC2_indexes, 1999 for PC1_indexes, and 2011 and 2012 for PC2_indexes and PC1_indexes, respectively (Table 6). The CC showed discontinuity in the ecosystem indicators data in 1990, 1999, and 2012 (Table 6).
<figure><img src="image_15.png" /><figcaption>**FIGURE 7** | Time-series of modelled biomass data. Significant shifts detected by STARS analysis are displayed on the plots as solid and dashed vertical line for shifts significant at *p* ≤ 0.05 and ≤0.1, respectively. Abbreviations, used in the figure, are explained in **Supplementary Table 4**.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**FIGURE 8** | Time-series of ecosystem indicators data. Significant shifts detected by STARS analysis are displayed on the plots as solid and dashed vertical line for shifts significant at *p* ≤ 0.05 and ≤0.1, respectively. Abbreviations, used in the figure, are explained in **Supplementary Table 4**.</figcaption></figure>

## DISCUSSION

This study provides a comprehensive overview of the changes observed in the PCSE in the last three decades (Figure 9). Ecosystem changes were assessed by using a large dataset that included information on ecosystem drivers and ecosystem state. Observed changes were confirmed by the detection of significant shifts in all used datasets.

## Description of PCSE Changes
### Changes in Ecosystem Drivers

The major changes in the ecosystem drivers occurred in the mid-1990s. Firstly, we observed a decrease in the fishing mortality of most functional groups. This can be related to the substantial reduction of the Portuguese fleet that decreased by 30% between 1989 and 1996 (Instituto Nacional de Estatística, 1986-1996). Secondly, in the mid-90s, an increase in SST, AMO, and EA, has
<table><caption>**TABLE 6** | Results of the sequential *t*-test analysis of regime shift (STARS) on principal component analysis (PCA) output and chronological clustering (CC) analysis.</caption><thead><tr><th rowspan="2">Data group</th><th colspan="4">Analysis</th></tr><tr><th colspan="3">STARS on PCs</th><th>CC</th></tr><tr><th></th><th>PC1</th><th>PC2</th><th>PC3</th><th>RSI</th></tr></thead><tbody><tr><td rowspan="2">Fishing mortality (FM)</td><td>1999↓</td><td>1995↓</td><td></td><td>-0.9/-0.6</td><td>1992</td></tr><tr><td></td><td>2003↑</td><td>2003↓</td><td>-0.4/-0.6</td><td>2003</td></tr><tr><td rowspan="2">Environmental/climatic (ENV_CLIM)</td><td>2012↓</td><td>2016↓</td><td>2014↑</td><td>-0.7/-0.4/2.2</td><td>-</td></tr><tr><td>1997↑</td><td></td><td></td><td>1.0/</td><td>1995<br>2002</td></tr><tr><td rowspan="2">Observed biomass (B_ob)</td><td>2016↑</td><td>2014↑</td><td></td><td>0.7/0.7</td><td>2009</td></tr><tr><td>1993↓</td><td>-</td><td></td><td>-1.3</td><td>1990</td></tr><tr><td rowspan="2">Modelled biomass (B_model)</td><td>2005↓</td><td>2003↑</td><td>-</td><td>-2.9/0.8</td><td>-</td></tr><tr><td>2013↓</td><td>2015↑</td><td>-</td><td>-0.1/0.4</td><td>-</td></tr><tr><td rowspan="3">Ecosystem indicators (Indices)</td><td>2008↓</td><td>1992↑</td><td></td><td>1.0</td><td>1998</td></tr><tr><td>2015↓</td><td></td><td></td><td>-2.1</td><td>-</td></tr><tr><td></td><td></td><td></td><td>-1.6</td><td>2015</td></tr><tr><td rowspan="2">Ecosystem indicators (Indices)</td><td>1999↓</td><td>1992↓</td><td></td><td>0.6/-0.4</td><td>1990<br>1999</td></tr><tr><td>2012↑</td><td>2011↑</td><td></td><td>1.8/1.4</td><td>2012</td></tr></tbody></table>

STARS regime shift years were identified in time series of PC1, PC2 and PC3 scores using *p* ≤ 0.05 significance level. Regime shift index (RSI) values for PC1/PC2/PC3 are presented. Symbols ↑ and ↓ indicate increasing and decreasing trend, respectively.

been reported. Around 1995, AMO entered a warm phase and since then the yearly SST means in North Atlantic including, the Iberian upwelling region has been ≥0.5°C higher than in the previous 65 years (Alheit et al., 2014). Later in the mid-2010s, we detected the change in winter NAO (that enters its high phase) and in the intensity of the upwelling. Changes observed in both these variables are in line with the findings by Borges et al. (2003) that demonstrated that in the Iberian shelf, high NAO phase is associated with more favourable upwelling winter wind conditions.

## Changes in Ecosystem State (Biomass)

During the studied period, we have seen considerable changes in the biomass of many functional groups. In the mid-90s the biomass of sardine, the main commercial pelagic fish started to decline while an abundance of other SPFs such as bogue and chub mackerel increased. Later, in the mid-2000s further changes in the biomass of functional groups in the ecosystem have been reported. During this time, we observed not only a further decrease of sardine and increase of other pelagics but also an increase of lower TL groups (zooplankton and invertebrates), sparids and cephalopods started to be apparent. After 2010, we still observed a decrease in sardine with a concurrent increase of other SPF and lower TL groups. Additionally, during this time an increase in predatory, high TL and demersal fish has been detected.

## Changes in Ecosystem State (Indicators)

Changes in the ecosystem were also reflected in the ecosystem indicators. The major changes in ecosystem indicators data were observed in the mid-90s and later in the mid-2010s.

In the mid-90s a major change was described by a decrease in Dem/Pel B, and catch related indices (Tot C, Dem C, Dem/Pel C, and Pred C). A decrease in Dem/Pel B can be explained by an increase in biomass of pelagic species, discussed above. On the other hand, a decrease in catch related indices is probably related to the reduction of fishing effort observed in the PCSE during the last three decades (Bueno-Pardo et al., 2020). Later in the mid-2010s, a further reduction in Dem/Pel B was observed. At this time, an increase in TL indices indicating a rise of higher TL fish has also been documented (mTLc, mTLco, Pred B, and Tot P/Tot B). Moreover, around this time, the indices that describe ecosystem stability and resilience shifted (FCI and R).

## Causes of Ecosystem Changes

The correlation analysis indicates that the changes in ecosystem trends (biomass and indicators) are statistically correlated with the changes in fishing and environmental/climatic variables.

### Fishing

Fishing has been documented as a major driver of the marine ecosystems worldwide (e.g., Pauly et al., 1998; Coll et al., 2008; Llope et al., 2011). Moreover, it has been highlighted as an important driver of the PCSE (Szalaj et al., 2021).
<figure><img src="image_17.png" /></figure>

The important role of fishing detected in this study can be explained by the considerable changes observed in the Portuguese fishing sector in the last 50 years. Since the late 1980s, reported landings dropped approximately 52%, from 197,000 to 94,700 t in 2014 (Bueno-Pardo et al., 2020). Moreover, after the implementation of the European Common Fisheries Policy (EU CFP) (since 1983), a steady decrease of fishing effort started to be notable with a decrease of 75 and 50% in the number of boats and fishers, respectively (Leitão et al., 2014b). A decrease in fishing effort, observed after enforcement of EU CFP measures, could affect an increase in biomass of some species observed in the study. For example, an increase in predatory higher TL fish, that manifested after 2010 can be linked with the reduction in fishing pressure. Observed in this analysis, a decrease in fishing mortality of most functional groups since the mid-90s (demersal

and benthopelagic piscivorous fish, rays, sparids, bogue, horse
mackerel, and mackerel), might have positively impacted the
high TL predatory fish and demersal fish whose abundance
increased and shifted after 2015 (with a time lag). A link between
fishing and high TL predatory fish established in this study
is based not only on the biomass data (observed increase in
the biomass of demersal and benthopelagic piscivorous fish)
but also based on ecosystem indicators. An increase in indices,
mTLc and mTLco, that are linked with the decreasing fishing
pressure (Pauly et al., 1998; Rochet and Trenkel, 2003), were
observed in this study. This is in line with the study by
Leitão (2015a) that assessed the profiles of Portuguese fisheries
landings since 1950. In this study, Leitão (2015a) documented
an increase in the mean TL of the catch, by approximately 0.2
units per decade.
However, in this study, we used fishing mortality data for only 10 functional groups, and this considerably limits the assessment of the impact of fishing on this system. Fishing mortality data for eight commercial groups, that were included in the biomass data (chub mackerel, sparids, flatfish, demersal and benthopelagic invertivorous fish, cephalopods, and sharks) were not considered as fishing drivers.

## Environmental/Climatic Drivers

The PC trends extracted from environmental/climatic drivers data showed significant correlations with all PC trends extracted from ecosystem data.

The importance of environmental variability to drive marine ecosystems has been reported worldwide (e.g., Parsons and Lear, 2001; Lehodey et al., 2020). This has been particularly prominent in the large upwelling systems (e.g., Shannon et al., 1976). The PCSE is located at the northern limit of the Canary Current upwelling system and during spring and summer it is affected by the seasonal upwelling (Wooster et al., 1976; Fiúza et al., 1982). Moreover, changes in the PCSE have been linked to environmental variability (e.g., Szalaj et al., 2021).

A decrease in sardine biomass and an increase in other pelagic fish and rays were both represented by PC1 in both data sets (observed and modelled biomass data groups). Additionally, in the observed biomass dataset, PC1 was represented by sparids, benthopelagic invertivorous fish, and tunas (for which data was not introduced in PCA_B_model) while for modelled biomass data set, PC1 was also represented by cephalopods. PC1_B_obsm and PC1_B_model showed a statistically significant correlation with PC1_ENV_CLIM and PC2_ENV_CLIM, respectively. PC1_ENV_CLIM represented SST, AMO and EA variables, while PC2_ENV_CLIM represented NAO and UI trends. Therefore, based on this analysis it seemed that all these environmental/climatic variables might have played a role in the overall biomass changes observed in the study.

It is possible that the changes in environmental drivers that were apparent in the mid-90s (SST, AMO, EA) affected the ecosystem and manifested with different time lag depends on the functional groups. In the PCSE, the productivity of species representing the functional groups for which significant shifts were observed in the mid-2000s, has been already documented as being influenced by environmental variability. Firstly, we have a large number of studies that link sardine recruitment and productivity with SST, chlorophyll concentration (Chla) (Santos et al., 2012; Leitão et al., 2014a; Garrido et al., 2017), NAO, wind conditions (Santos et al., 2012; Leitão et al., 2014a) and coastal upwelling (Santos et al., 2001; Santos et al., 2012; Malta et al., 2016).

Moreover, other pelagic species such as Atlantic horse mackerel (*T. trachurus*) and anchovy (*Engraulis encrasicolus*) catch rates were linked with upwelling, the magnitude of northern wind and SST (Santos et al., 2001; Ullah et al., 2012; Leitão, 2015b). Also, among cephalopods, cuttlefish (*Sepia oficinalis*) catch rates were associated with the SST, winter NAO and autumn wind magnitude (Ullah et al., 2012). Furthermore, by using a statistical multi-model approach, Baptista et al. (2015) found that environmental factors such as winter SST and NAO

affected short term variation of catch rates of common two-banded seabream (*Diplodus vulgaris*) in the south and southwest of Portugal. Also, they found that these species catch rates were associated with autumn and a yearly wind magnitude component, that can be treated as a proxy for the upwelling condition. All these studies were performed using statistical models. It is therefore encouraging that results from the study using different methods such as the outputs of ecosystem model and discontinuity analyses, established similar associations.

After 2010, we still observed a decrease in sardine and a concurrent increase in other SPF and lower TL groups. Additionally, during this time an increase in predatory, higher TL and demersal fish was also detected. When analysing shifts in data describing drivers, we observe that around this time the shift in winter NAO (that enters its high phase) and in the intensity of upwelling was detected. In the Iberian shelf, a high NAO phase is associated with more favourable upwelling winter wind conditions (Borges et al., 2003) that may favour lower trophic groups (phytoplankton, zooplankton, benthos). This can quickly promote the abundance of lower TL groups that rely directly on plankton as a food source (Barth et al., 2007). An increase in lower TL groups was initiated earlier in the mid-2000s but the shift in NAO in the mid-2010s might have further enhanced the change that initiated earlier.

## Internal Drivers (Trophic Interactions)

The effect of environmental variability on some species (functional groups) might have been intensified by the internal, trophic interactions that exist between these species. The importance of trophic interactions in the dynamics of the marine ecosystems has been well documented worldwide (e.g., Coll et al., 2008; Alexander et al., 2015). Also, it has been demonstrated as the main factor driving PCSE dynamics between 1986 and 2017 (Szalaj et al., 2021).

In the context of this study, an increase in lower TL groups due to the environmental conditions, could pronounce trough the food web and promote the abundance of fish of prey that consequently might have benefited the groups at higher TL including cephalopods, demersal fish, and higher TL predatory fish. Also, the effect of environmental variability on some species (functional groups) might be exaggerated by the existence of the predator-prey loop between these species. This can be a case for sardine and chub mackerel. The increase in the biomass of chub mackerel observed after 2001 (that possibly benefited from increased temperature) may have caused increased pressure on sardine recruitment as it has been established that sardine eggs are chub mackerel's preferred prey (Garrido et al., 2015). At the same time, increased biomass of other species such as bogue, horse mackerel and mackerel might have added pressure on sardine due to the predation on its eggs and juveniles (Garrido et al., 2015).

## Implications to Ecosystem Structure and Function (ENA Indices)

Other than the changes evident from the direct observations and modelled biomass, the ENA indices provide information about the changes that occurred at the whole ecosystem level.
Analysis of ENA indices showed that not only has the abundance of functional groups altered but also the cycles and pathways that are implied by the ecosystem food web showed considerable alterations within the last thirty years.

ENA indices (R, FCI, and A/C), indicate that after 2010, the ecosystem became less stable and less resilient with reduced maturity and complexity. Lower stability and resilience, reflected by the R, is implied by the reduction in the number of pathways through which energy flow is transferred. Moreover, reduction in R and increase in A/C implies lower maturity of the system which become more constrained and less complex. Lower FCI indicates that the system's recycling and its ability to conserve nutrients are reduced and therefore, its capacity to endure in the face of resource scarcity has diminished. These characteristics indicate structural differences (Finn, 1976) and decreased stability (DeAngelis et al., 1989; Loreau, 1994). Reduced resilience and stability concur with the observed decrease in diversity indices as it has been demonstrated that a system with greater biological diversity is usually more resilient than the one with less diversity (Oliver et al., 2015). It is not clear, what has caused the changes in the ecosystem structure and its flows. However, the most prominent change (the change of the highest magnitude) observed in the last decades was the decline of sardine. Indeed, the effects of the changes in dominant SPF on the ecosystem food web in the upwelling systems, have been broadly discussed (Cury et al., 2000). Moreover, based on the Ecopath model and description of the Iberian shelf food web structure, it has been demonstrated that sardine played an important role in the Iberian ecosystem, being an important key linkage between primary producers and top predators (Veiga-Malta et al., 2019). Moreover, the ecosim simulations demonstrated that sardine has an important role in a wasp-waist scenario (Szalaj et al., 2021). This confirmed its key role in the Iberian upwelling system as a biomass-dominant species at mid-TLs that channel most of the energy flow in the ecosystem. Therefore, it seems probable that sardine decline disturbed important trophic links and pathways through which the flow of energy passes. At the same time, the other SPF were not well linked to the rest of the food web with few animals predating on them and being not very important for the fishery. For example, bogue was only caught as a bycatch of the purse seine fishery and catches of chub mackerel only become important after 2000 when a decline of sardine occurred and fisheries managed to find alternatives for economic compensation (Instituto Nacional de Estatística, 2000-2018). Therefore, when pelagic energy flow was redirected to the other SPF (that were not as important prey for top predators and fisheries as sardine was before 1995), disruption in the energy flow pathways occurred. Consequently, ecosystem stability was reduced, and the ecosystem became less resilient. Similar changes were described in the Northern Benguela upwelling system, where a shift in dominant pelagic species coincided with the simplification of the food web and reduced resilience (Heymans and Tomczak, 2016).

## Assessment of the Methodology

Methods that comprised of ITA methodology has been used to detect ecological and climatic regime shifts and they have

been widely discussed and considered adequate (e.g., Rodionov and Overland, 2005; Möllmann et al., 2009; Diekmann and Möllmann, 2010; Tomczak et al., 2013). However, the use of EwE outputs in this type of analysis is rare and to our knowledge there are only two studies that applied this approach to detect and describe regime shifts. One in the Baltic Sea and the other in the Northern Benguela (Tomczak et al., 2013; Heymans and Tomczak, 2016). The main advantage of this approach is the possibility to use more data inputs, that EwE provides, for example, the data of functional groups for which observations do not exist and ecosystem indicators. Among the latter, ENA indicators can be especially valuable in the studies of ecosystem regime shifts as discussed by Tomczak et al. (2013). For example, the redundancy that was linked with ecosystem resilience (Heymans et al., 2007) is directly related to regime shift definition. Therefore, we argue that using this approach is very beneficial in the studies that aim to detect and describe ecosystem regime shifts.

However, we are aware of the limitations of our analyses, such as high cross and auto-correlations. Planque and Arneberg (2017) criticised the use of PCA as a method for ecosystem assessment of marine ecosystems because the ITAs inputs usually share a similar trend and have autocorrelation properties. The authors demonstrated that most of the patterns revealed by the PCA can emerge from random time-series and that the fraction of the variance that cannot be accounted by random processes is minimal. In the context of our study these findings are worrying and this limitation should be taken into consideration when interpreting the results. We tried to minimalize this effect by excluding highly correlated variables prior to PCA. Also, the ITA methodology used in this study combines different methods and not all of them rely on PCA. For example, STARS performed on separate time series accounts for autocorrelation (prewhitening method described in the "Material and Methods" section) and detected shifts in similar years (time ranges) that STARS performed on PCA output. Moreover, the CC method confirmed dissimilarity in the datasets in similar years as STARS on PCA results. However, to confirm and consolidate the conclusions of this study, as a future recommendation we would suggest to replicate the analysis using other methods that are not so sensitive to cross and autocorrelation and compare its outputs with the results obtained in this study.

Moreover, our analysis is fed with data that constitute the output of a modelling approach therefore our results are affected by the assumptions of the model. EwE is a commonly used approach (Colléter et al., 2015) and its best practices, pros and cons have been widely discussed (Heymans et al., 2016; Walters et al., 1997). When building and calibrating the EwE model for the PCSE, we considered the best practice guide (Heymans et al., 2016), reported uncertainties and thoroughly discussed limitations (Szalaj et al., 2021). Moreover, the model was well fitted to most groups with observation data, and we have confidence in the model and data it provided (Veiga-Malta et al., 2019; Szalaj et al., 2021). This was also reflected in this study because the results obtained from the analyses performed on observed biomass are in line with the results obtained from modelled biomass.
## Regime Shifts in the PCSE and North Atlantic

The performed analyses allowed to detect three shifts in data, in the mid-90s, mid-2000s, and mid-2010s. However, based on the analysis of the results and the review of the literature, only two of them, in the mid-90s and mid-2010s, could be considered as shifts to a new regime.

The shifts in the mid-90s and mid-2010s were detected by STARS on all PC trends regardless of the data group used, while the shift in the mid-2000s was mainly detected in the biomass data (when considering main PCs). Moreover, the timing of the shifts detected in the mid-90s and mid-2010s in this study coincided with the time of ecosystem regime shifts documented in the other regions of the North Atlantic. The mid-90s shift is well documented in the northern hemisphere and a large set of previous studies described it and linked it with large scale atmospheric variations (e.g., Cabrero et al., 2018; Alheit et al., 2014; Beaugrand et al., 2015). For example, Alheit et al. (2014), linked the complex ocean-atmosphere changes in the mid-90s with a regime shift in the ecosystems of the Northeastern and Central Atlantic. He discussed that when in the mid-1990s, three events including (1) decrease of the NAO index, (2) the weakening of the Atlantic Meridional Overturning Circulation (AMOC), and (3) the contraction of the subpolar gyre (SPG) occurred, a substantial redistribution of water masses in the Northeast Atlantic happened. This affected most ecosystems in the region causing a change in the abundance of the small pelagic clupeoid fish and their northward migration. Moreover, Cabrero et al. (2018) demonstrated that in the Iberian shelf, sardine entered a new regime in the late 1990s triggered by large-scale ocean climate processes. He concluded that the combination of climate and oceanographic conditions, including NAO, EA, SST, AMO, and UI intensity, contributed to this regime shift.

Also, a reference to the mid-2010s shift in the Northern Atlantic appeared in the recent literature. For example, Tanner et al. (2020) found that the abrupt changes in the growth of two deep-sea scorpaenid fishes off the central islands of the Azores and the island of Madeira occurred not only in the late 90s but also around 2014. This coincides with a significant shift in the mid-2010s detected in this study. Tanner et al. (2020) associated the shift in 2014 with a large-scale abrupt biological shift that was predicted by Beaugrand et al. (2019) to occur after 2014 in response to a strong El Nino event that affected the Northern Hemisphere climate. In the study by Beaugrand et al. (2019), a model based on the numerically generated thermal niche of pseudo-communities was used and the shifts predicted by the model were validated with observed community shifts in 14 oceanic regions. This study is of special interest to us, as the shift in the mid-2010s, that was detected in our analysis, is not well documented in the scientific literature. A few substantial climate events occurred around this time: exceptionally cold North Atlantic Ocean temperatures, 2015 severe European heatwave (Duchez et al., 2016) and the change in AMO phase (Frajka-Williams et al., 2017). However, these events are quite recent, and their ecological consequences are not yet well studied and understood.

The mid-2000 regime shift that was detected in this study did not relate to any relevant environmental pressure observed in the North Atlantic in our dataset or the literature. Generally, the time between 2005 and 2014 is considered a relatively quiet period when referring to the climatic divers in the Northern hemisphere (Beaugrand et al., 2015). Moreover, this shift was not reflected in the changes observed in the ecosystem indicators analysed in the study. Therefore, we considered the period between the first shift in the mid-90s and the shift in the mid-2010s as the “transition regime” during which the changes triggered by a previous shift in the mid-90s manifested in the ecosystem with the time lag whose length depends on the functional group.

From this analysis it is not clear whether the shift was caused solely by the environmental/climatic causes because both analysed ecosystem drivers' trends (fishing and environmental/climatic) correlated with ecosystem state trends. However, based on the literature and a well-documented climatic shift in the North Atlantic we can infer that climatic variability played a major role. However, the effect of fishing pressure should not be completely disregarded as fishing is a well-documented driver of marine ecosystems worldwide (e.g., Coll et al., 2008; Alexander et al., 2015) and substantial changes in the fishing sector has been documented during the last three decades in the PCSE (Leitão et al., 2014b; Bueno-Pardo et al., 2020). Moreover, regime shifts are commonly triggered by a combination of fishing and environmental/climatic factors acting together (e.g., Daskalov et al., 2007; Möllmann et al., 2009). In the case of the PCSE, although we have documented a decrease in fishing effort (Bueno-Pardo et al., 2020), many Portuguese fish stocks are considered overfished (Leitão, 2015a). As we discussed in section “Fishing,” the impact of fishing on the ecosystem has not been fully captured in this study. Therefore, the hypothesis of the fishing cause of the PCSE regime shift should be explored in future studies.

At this stage, it is impossible to say if the current regime is stable or reversible as the PCSE as an upwelling system might have a wide range of natural variability (Black et al., 2014). If the observed ecosystem food web changes are linked with sardine decline as discussed above (in the previous subsection) then it is unknown what effect sardine recovery may have on the ecosystem. The latest ICES sardine stock assessment indicated that the 2019 sardine recruitment estimate constitutes the highest value since 2004 and is above the long-term geometric mean (ICES, 2020). That might suggest that the sardine productivity trend is shifting. However, at that stage it is unknown whether the increase in sardine recruitment will continue and sardine stock will recover to the level prior to the mid-1990s and what impact would that potential increase have on the food web.

## Implications for Management

The potential occurrence of regime shifts might have considerable implications for the productivity of the commercial stocks and consequently to the fisheries. Therefore, the awareness of regime shift in the management of the ecosystems is of utmost importance and should be included in the management plan of marine systems (deYoung et al., 2008). This management should include not only timely detection and adaptation to a new
regime but also prevention (deYoung et al., 2008; Crépin et al., 2012). The latter is crucial and should focus on building and maintaining the resilience of desired ecosystem states (Levin and Lubchenco, 2008). Nevertheless, sometimes, regime shifts cannot be prevented, and in these circumstances, timely detection is important in order to act quickly to alleviate the adverse effects and adapt to a new reality (Crépin et al., 2012).

After the detection of regime shifts, the adaptation in the form of adaptive management should take place. For example, a regime-specific harvest rate strategy could be implemented when an indication toward a shift to an alternate ecosystem state exists (Polovina, 2005; MacCall, 2002). This study demonstrated that one of the possible mechanism behind the described regime shift is environmentally related therefore in the context of PCSE management, efforts should be focused on establishing reference points at which the change in environmental pressure might result in an abrupt change in the ecosystem (Large et al., 2015).

## CONCLUSION

The ITA presented in this study pointed to the occurrence of two regime shifts in the PCSE in the mid-90s and around the mid-2010s. Interpretation of results allowed for the identification of three periods: “early regime” till the mid-90s, “transition regime” till mid-2010s, and “late regime” since the mid-2010s in the last three decades. Differences between “early regime” and “late regime” are evident in all components of the ecosystem and at multiple TLs. They include changes in species abundance as well as community composition. Moreover, ENA indices indicate that the internal food web structure and organisation of the ecosystem have changed. Specifically, detected regime shifts are characterised by changes in the pelagic community: (1) it became more dominant when compared to the demersal community; (2) a shift from the main commercial stock of sardine to less commercially viable small pelagics had occurred. Moreover, an increase in higher TL predatory fish and demersal fish was observed.

Changes described in the ecosystem were probably caused by two types of external drivers. First, environmental forcing (related to an increase in SST, AMO, EA, NAO, and intensification of the upwelling conditions) possibly affected lower TL groups and SPF, and second, reduction in fishing pressure probably benefited high

TL predatory and demersal fish. Moreover, the internal triggers, related to indirect trophic interactions also might play a role in the observed ecosystem change. An increased abundance of lower TL functional groups and fish of prey might benefit higher TL groups such as cephalopods, demersal and predatory fish. Furthermore, the competition and intraguild predation between sardine and other pelagic fish might have contributed to the changes in the pelagic community.

## DATA AVAILABILITY STATEMENT

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

## AUTHOR CONTRIBUTIONS

DS conceptualized the study, wrote the original draft, conducted the modelling, performed statistical analyses, and interpreted the results. AS, PR, and HC conceptualized the study, contributed to the interpretation of results, reviewed the manuscript, and assisted in funding acquisition. All authors contributed to the article and approved the submitted version.

## FUNDING

This publication was financed by Portuguese national funds through FCT – Fundação IP under project reference UIDB/04292/2020, and by the European Union’s Horizon 2020 Research and Innovation Programme under grant agreement N810139: Project Portugal Twinning for Innovation and Excellence in Marine Science and Earth Observation – PORTWIMS. Work of DS was funded by FCT (Fundação para a Ciência e a Tecnologia) Ph.D. fellowship (ref PD/BD/114540/2016).

## SUPPLEMENTARY MATERIAL

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

## REFERENCES

Alexander, K. A., Heymans, J. J., Magill, S., Tomczak, M. T., Holmes, S. J., and Wilding, T. A. (2015). Investigating the recent decline in gadoid stocks in the west of Scotland shelf ecosystem using a foodweb model. *ICES J. Mar. Sci.* 72, 436–449. doi: 10.1093/icesjms/fsu149

Alheit, J., Licandro, P., Coombs, S., Garcia, A., Giráldez, A., Santamaría, M. T. G., et al. (2014). Atlantic Multidecadal Oscillation (AMO) modulates dynamics of small pelagic fishes and ecosystem regime shifts in the eastern North and Central Atlantic. *J. Mar. Syst.* 133, 88–102. doi: 10.1016/j.jimarsys.2013.11.002

Baptista, V., Campos, C. J. A., and Leitão, F. (2015). The influence of environmental factors and fishing pressure on catch rates of *Diplodus vulgaris*. *Estuaries Coasts* 39, 258–272. doi: 10.1007/s12237-015-9990-y

Baptista, V., Silva, P. L., Relvas, P., Teodósio, M. A., and Leitão, F. (2017). Sea surface temperature variability along the Portuguese coast since 1950. *Int. J. Climatol.* 38, 1145–1160. doi: 10.1002/joc.5231

Barth, J. A., Menge, B. A., Lubchenco, J., Chan, F., Bane, J. M., Kirincich, A. R., et al. (2007). Delayed upwelling alters nearshore coastal ocean ecosystems in the northern California current. *Proc. Natl. Acad. Sci. U.S.A.* 104, 3719–3724. doi: 10.1073/Proc Natl. Acad. Sci. U.S.A. 0700462104

Barton, E. D., Field, D. B., and Roy, C. (2013). Canary current upwelling: more or less? *Prog. Oceanogr.* 116, 167–178. doi: 10.1016/j.pocean.2013.07.007

Beaugrand, G., Conversi, A., Atkinson, A., Cloern, J., Chiba, S., Fonda-Umani, S., et al. (2019). Prediction of unprecedented biological shifts in the global ocean. *Nat. Clim. Change* 9, 237–243. doi: 10.1038/s41558-019-0420-1

Beaugrand, G., Conversi, A., Chiba, S., Edwards, M., Fonda-Umani, S., Greene, C., et al. (2015). Synchronous marine pelagic regime shifts in the Northern Hemisphere. *Philos. Trans. R. Soc. B* 370:20130272. doi: 10.1098/rstb.2013.0272

Black, B. A., Sydeman, W. J., Frank, D. C., Griffin, D., Stahle, D. W., Garcia-Reyes, M., et al. (2014). Six centuries of variability and extremes in a coupled marine-terrestrial ecosystem. *Science* 345, 1498–1502. doi: 10.1126/science.1253209
Borges, M. F., Santos, A. M. P., Crato, N., Mendes, H., and Mota, B. (2003). Sardine regime shifts off Portugal: a time series analysis of catches and wind conditions. *Sci. Mar.* 67, 235-244. doi: 10.3989/scimar.2003.67s1235

Bueno-Pardo, J., Pierce, G. J., Cabecinha, E., Grilo, C., Assis, J., Valavanis, V., et al. (2020). Trends and drivers of marine fish landings in Portugal since its entrance in the European Union. *ICES J. Mar. Sci.* 77, 988-1001. doi: 10.1093/icesjms/fsaa010

Cabrero, Á, González-Nuevo, G., Gago, J., and Cabanas, J. M. (2018). Study of sardine (*Sardina pilchardus*) regime shifts in the Iberian Atlantic shelf waters. *Fish. Oceanogr.* 28, 1-12. doi: 10.1111/fog.12410

Christensen, V. (1995). Ecosystem maturity – towards quantification. *Ecol. Model.* 77, 3-32. doi: 10.1016/0304-3800(93)E0073-C

Christensen, V., and Pauly, D. (1992). ECOPATH II – a software for balancing steady-state ecosystem models and calculating network characteristics. *Ecol. Model.* 61, 169-185. doi: 10.1016/0304-3800(92)90016-8

Christensen, V., and Walters, C. J. (2004). Ecopath with Ecosim: methods, capabilities and limitations. *Ecol. Model.* 172, 109-139. doi: 10.1016/j.ecolmodel.2003.09.003

Coll, M., Palomera, I., Tudela, S., and Dowd, M. (2008). Food-web dynamics in the South Catalan Sea ecosystem (NW Mediterranean) for 1978–2003. *Ecol. Model.* 217, 95-116. doi: 10.1016/j.ecolmodel.2008.06.013

Colléter, M., Valls, A., Guitton, J., Gascuel, D., Pauly, D., and Christensen, V. (2015). Global overview of the applications of the Ecopath with Ecosim modeling approach using the EcoBase models repository. *Ecol. Model.* 302, 42-53. doi: 10.1016/j.ecolmodel.2015.01.025

Collie, J. S., Richardson, K., and Steele, J. H. (2004). Regime shifts: can ecological theory illuminate the mechanisms? *Prog. Oceanogr.* 60, 281-302. doi: 10.1016/j.pocean.2004.02.013

Crépin, A. S., Biggs, R., Polasky, S., Troell, M., and de Zeeuw, A. (2012). Regime shifts and management. *Ecol. Econ.* 84, 15-22. doi: 10.1016/j.ecolecon.2012.09.003

Cury, P. M., Shannon, L. J., Roux, J. P., Daskalov, G. M., Jarre, A., Moloney, C. L., et al. (2005). Trophodynamic indicators for an ecosystem approach to fisheries. *ICES J. Mar. Sci.* 62, 430-442. doi: 10.1016/j.icesjms.2004.12.006

Cury, P., and Shannon, L. J. (2004). Regime shifts in upwelling ecosystems: observed changes and possible mechanisms in the northern and southern Benguela. *Prog. Oceanogr.* 60, 223-243. doi: 10.1016/j.pocean.2004.02.007

Cury, P., Bakun, A., Crawford, R. J. M., Jarre, A., Quiñones, R. A., Shannon, L. J., et al. (2000). Small pelagics in upwelling systems: Patterns of interaction and structural changes in “wasp-waist” ecosystems. *ICES J. Mar. Sci.* 57, 603-618. doi: 10.1006/jmsc.2000.0712

Daskalov, G. (2002). Overfishing drives a trophic cascade in the Black Sea. *Mar. Ecol. Prog. Ser.* 225, 53-63. doi: 10.3354/meps225053

Daskalov, G. M., Grishin, A., Rodionov, S., and Mihneva, V. (2007). Trophic cascades triggered by overfishing reveal possible mechanisms of ecosystem regime shifts. *Proc. Natl. Acad. Sci. U.S.A.* 104, 10518-10523.

DeAngelis, D., Mulholland, P., Palumbo, A., Steinman, A., Huston, M., and Elwood, J. (1989). Nutrient dynamics and food-web stability. *Annu. Rev. Ecol. Syst.* 20, 71-95. doi: 10.1146/annurev.es.20.110189.000443

deYoung, B., Barange, M., Beaugrand, G., Harris, R., Perry, I. R., Scheffer, M., et al. (2008). Regime shifts in marine ecosystems: detection, prediction and management. *Trends Ecol. Evol.* 23, 402-409. doi: 10.1016/j.tree.2008.03.008

deYoung, B., Harris, R., Alheit, J., Beaugrand, G., Mantua, N. J., and Shannon, L. J. (2004). Detecting regime shifts in the ocean: data considerations. *Prog. Oceanogr.* 60, 143-164. doi: 10.1016/j.pocean.2004.02.017

Diekmann, R., and Möllmann, C. (2010). Integrated ecosystem assessment of seven Baltic Sea areas covering the last three decades. *ICES Cooperat. Res. Rep.* 302:90. doi: 10.17895/ices.pub.5413

Diekmann, R., Otto, S. A., and Möllmann, C. (2012). "Towards integrated ecosystem assessments (IEAs) of the Baltic Sea: investigating ecosystem state and historical development," in *Climate Impacts on the Baltic Sea: From Science to Policy*, eds M. Reckermann, K. Brander, B. R. MacKenzie, and A. Omstedt (Heidelberg: Springer Earth System Sciences Series), 161-201.

Duchez, A., Frajka-Williams, E., Josey, S. A., Evans, D. G., Grist, J. P., Marsh, R., et al. (2016). Drivers of exceptionally cold North Atlantic Ocean temperatures and their link to the 2015 European heat wave. *Environ. Res. Lett.* 11:074004. doi: 10.1088/1748-9326/11/7/074004

Finn, J. T. (1976). Measures of ecosystem structure and function derived from analysis of flows. *J. Theor. Biol.* 56, 363-380. doi: 10.1016/S0022-5193(76) 80080-X

Fiúza, A. F. D., de Macedo, M. E., and Guerreiro, M. R. (1982). Climatological space and time variation of the Portuguese coastal upwelling. *Oceanol. Acta* 5, 31-40.

Folke, C., Carpenter, S., Walker, B., Scheffer, M., Elmqvist, T., Gunderson, L., et al. (2004). Regime shifts, resilience, and biodiversity in ecosystem management. *Annu. Rev. Ecol. Evol. Syst.* 35, 557-581. doi: 10.1146/annurev.ecolsys.35.021103.105711

Frajka-Williams, E., Beaulieu, C., and Duchez, A. (2017). Emerging negative Atlantic Multidecadal Oscillation index in spite of warm subtropics. *Sci. Rep.* 7:11224. doi: 10.1038/s41598-017-11046-x

Garrido, S., Silva, A., Marques, V., Figueiredo, I., Bryère, P., Mangin, A., et al. (2017). Temperature and food-mediated variability of European Atlantic sardine recruitment. *Prog. Oceanogr.* 159, 267-275. doi: 10.1016/j.pocean.2017.10.006

Garrido, S., Silva, A., Pastor, J., Dominguez, R., Silva, A. V., and Santos, A. M. (2015). Trophic ecology of pelagic fish species off the Iberian coast: Diet overlap, cannibalism and intraguild predation. *Mar. Ecol. Prog. Ser.* 539, 271-286. doi: 10.3354/meps11506

Heymans, J. J., Coll, M., Link, J. S., Mackinson, S., Steenbeek, J., Walters, C., et al. (2016). Best practice in Ecopath with Ecosim food-web models for ecosystem-based management. *Ecol. Model.* 331, 173-184. doi: 10.1016/j.ecolmodel.2015.12.007

Heymans, J. J., Guenette, S., and Christensen, V. (2007). Evaluating network analysis indicators of ecosystem status in the Gulf of Alaska. *Ecosystems* 10, 488-502. doi: 10.1007/s10021-007-9034-y

Heymans, J. J., and Tomczak, M. T. (2016). Regime shifts in the Northern Benguela ecosystem: challenges for management. *Ecol. Model.* 331, 151-159. doi: 10.1016/j.ecolmodel.2015.10.027

Holling, C. S. (1973). Resilience and stability of ecological systems. *Annu. Rev. Ecol. Syst.* 4, 1-23. doi: 10.1146/annurev.es.04.110173.000245

ICES (2008). *ICES IBTSWG REPORT 2012 Report of the International Bottom Trawl Survey Working Group (IBTSWG)*. Lisbon: ICES, 27-30.

ICES (2016). *Second Interim report of the Working Group on Acoustic and Egg Surveys for sardine and anchovy in ICES Areas VII, VIII and IX (WGACEGG), 16-26 November 2015*. Lowestoft: ICES.

ICES (2017a). *Report of the Working Group on Southern Horse Mackerel, Anchovy and Sardine (WGHANSA), 24-29 June 2017*. Bilbao: ICES.

ICES (2017b). *Manual of the IBTS North Eastern Atlantic surveys. Series of ICES Survey Protocols SISP*, Vol. 15. Burnaby, BC: ICES, 92.

ICES (2018a). *Report of the Working Group for the Bay of Biscay and the Iberian Waters Ecoregion (WGBIE), 3-10 May 2018, ICES HQ*. Copenhagen: ICES.

ICES (2018b). *Report of the Working Group on Southern Horse Mackerel, Anchovy and Sardine (WGHANSA), 26-30 June 2018*. Lisbon: ICES.

ICES (2020). *Working Group on Southern Horse Mackerel, Anchovy and Sardine (WGHANSA). Draft report. ICES Scientific Reports*. Burnaby, BC: ICES, doi: 10.17895/ices.pub.5977

Instituto Nacional de Estatística (1986-1996). *INE - Instituto Nacional de Estatística, 1987-1997. Estatísticas da Pesca 2001-2018*. Lisboa: INE.

Instituto Nacional de Estatística (2000-2018). *INE - Instituto Nacional de Estatística, 2001-2017. Estatísticas da Pesca 2001-2018*. Lisboa: INE.

Juggins, S. (2020). *riocha: Analysis of Quaternary Science Data, R package version (0.9-26).*

Kempton, R., and Taylor, L. (1976). Models and statistics for species diversity. *Nature* 262, 818-820. doi: 10.1038/262818a0

Kendall, M. G. (1954). Note on bias in the estimation of autocorrelation. *Biometrika* 41, 403-404.

Kimura, D. K., Balsiger, J. W., and Ito, D. H. (1984). Generalized stock reduction analysis. *Can. J. Fish. Aquat. Sci.* 41, 1325-1333. doi: 10.1139/f84-162

Large, I. S., Fay, G., Friedland, K. D., and Link, J. S. (2015). Critical points in ecosystem responses to fishing and environmental pressures. *Mar. Ecol. Prog. Ser.* 521, 1-17. doi: 10.3354/meps11165

Lees, K., Pitois, S., Scott, C., Frid, C., and Mackinson, S. (2006). Characterizing regime shifts in the marine environment. *Fish Fish.* 7, 104-127. doi: 10.1111/j. 1467-2979.2006.00215.x
Legendre, P., Dallot, S., and Legendre, L. (1985). Succession of species within a community: chronological clustering, with applications to marine and freshwater zooplankton. *Am. Nat.* 125, 257-288.

Lehodey, P., Bertrand, A., Hobday, A. J., Kiyofuji, H., McClatchie, S., Menkès, C. E., et al. (2020). "ENSO impact on marine fisheries and ecosystems," in *El Niño Southern Oscillation in a Changing Climate*, eds M. J. McPhaden, A. Santoso, and W. Cai (Hoboken, NJ: John Wiley & Sons, Inc), 429-451. doi: 10.1002/9781119548164.ch19

Leitão, F. (2015a). Landing profiles of Portuguese fisheries: assessing the state of stocks. *Fish. Manag. Ecol.* 22, 152-163. doi: 10.1111/fme.12112

Leitão, F. (2015b). Time series analyses reveal environmental and fisheries controls on Atlantic horse mackerel (*Trachurus trachurus*) catch rates. *Continent. Shelf Res.* 111, 342-352. doi: 10.1016/j.csr.2015.08.026

Leitão, F., Alms, V., and Erzini, K. (2014a). A multi-model approach to evaluate the role of environmental variability and fishing pressure in sardine fisheries. *J. Mar. Syst.* 139, 128-138. doi: 10.1016/j.jmarsys.2014.05.013

Leitão, F., Baptista, V., Vieira, V., Laginha Silva, P., Relvas, P., and Teodósis, M. A. (2019). A 60-year time series analyses of the upwelling along the portuguese coast. *Water* 11:1285. doi: 10.3390/w11061285

Leitão, F., Baptista, V., Zeller, D., and Erzini, K. (2014b). Reconstructed catches and trends for mainland Portugal fisheries between 1938 and 2009: Implications for sustainability, domestic fish supply and imports. *Fish. Res.* 155, 33-50. doi: 10.1016/j.fishres.2014.02.012

Levin, P. S., and Möllmann, C. (2015). Marine ecosystem regime shifts: challenges and opportunities for ecosystem-based management. *Phil. Trans. R. Soc. B* 370:20130275. doi: 10.1098/rstb.2013.0275

Levin, S. A., and Lubchenco, J. (2008). Resilience, robustness, and marine ecosystem-based management. *Bioscience* 58, 27-32. doi: 10.1641/b580107

Llope, M., Daskalov, G. M., Rouyer, T. A., Mihneva, V., Chan, K.-S., Grishin, A. N., et al. (2011). Overfishing of top predators eroded the resilience of the Black Sea system regardless of the climate and anthropogenic conditions. *Glob. Chang. Biol.* 17, 1251-1265. doi: 10.1111/j.1365-2486.2010.02331.x

Loreau, M. (1994). Material cycling and the stability of ecosystems. *Am. Nat.* 143, 508-513.

MacCall, A. D. (2002). Fishery management and stock rebuilding prospects under conditions of low frequency variability and species interactions. *Bull. Mar. Sci.* 70, 613-628.

Mageau, M. T., Costanza, R., and Ulanowicz, R. E. (1998). Quantifying the trends expected in developing ecosystems. *Ecol. Model.* 112, 1-22. doi: 10.1016/S0304-3800(98)00092-1

Malta, T., Santos, P. T., Santos, A. M. P., Rufino, M., and Silva, A. (2016). Long-term variations in Ibero-Atlantic sardine (*Sardina pilchardus*) population dynamics: relation to environmental conditions and exploitation history. *Fish. Res.* 179, 47-56. doi: 10.1016/j.fishres.2016.02.009

Marriot, F. H. C., and Pope, J. A. (1954). Bias in the estimation of autocorrelations. *Biometrika* 41, 390-402.

Martins, M. M., Skagen, D., Marques, V., Zwolinski, J., and Silva, A. (2013). Changes in the abundance and spatial distribution of the Atlantic chub mackerel (*Scomber colias*) in the pelagic ecosystem and fisheries off Portugal. *Sci. Mar.* 77, 551-563. doi: 10.3989/scimar.03861.07B

May, R. M. (1977). Thresholds and breakpoints in ecosystems with a multiplicity of stable states. *Nature* 269, 471-477. doi: 10.1038/269471a0

McKinnell, S. M., Brodeur, R. D., Hanawa, K., Hollowed, A. B., Polovina, J. J., and Zhang, C. I. (2001). An introduction to the beyond El Nino conference: climate variability and marine ecosystem impact from the tropics to the Arctic. *Prog. Oceanogr.* 49, 1-6. doi: 10.1016/S0079-6611(01)00012-X

Miranda, P. M. A., Alves, J. M. R., and Serra, N. (2012). Climate change and upwelling: response of Iberian upwelling to atmospheric forcing in a regional climate scenario. *Clim. Dyn.* 40, 2813-2824. doi: 10.1007/s00382-012-1442-9

Möllmann, C., Diekmann, R., Muller-Karulis, B., Kornilovs, G., Plikshs, M., and Axe, P. (2009). Reorganization of a large marine ecosystem due to atmospheric and anthropogenic pressure: a discontinuous regime shift in the Central Baltic Sea. *Glob. Chang. Biol.* 15, 1377-1393. doi: 10.1111/j.1365-2486.2008.01184.x

Monaco, M. E., and Ulanowicz, R. E. (1997). Comparative ecosystem trophic structure of tree US mid Atlantic estuaries. *Mar. Ecol. Prog. Ser.* 161, 239-254.

Oksanen, J., Guillaume Blanchet, F., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., et al. (2019). *vegan: Community Ecology Package*. R package version 2.5-6.

Oliver, T. H., Heard, M. S., Isaac, N. J. B., Roy, B. D., Procter, D., Eigenbrod, F., et al. (2015). Biodiversity and resilience of ecosystem functions. *Trends Ecol. Evol.* 30, 673-684. doi: 10.1016/j.tree.2015.08.009

Parsons, L. S., and Lear, W. H. (2001). Climate variability and marine ecosystem impacts: a North Atlantic perspective. *Prog. Oceanogr.* 49, 167-188. doi: 10.1016/S0079-6611(01)00021-0

Pauly, D., Christensen, V., Dalsgaard, J., Froese, R., and Torres, F. (1998). Fishing down marine food webs. *Science* 279, 860-863. doi: 10.1126/science.279.5352.860

Pauly, D., Christensen, V., and Walters, C. J. (2000). Ecopath, Ecosim and Ecospace as tools for evaluating ecosystem impact of fisheries. *ICES J. Mar. Sci.* 57, 697-706. doi: 10.1006/jmsc.2000.0726

Pires, A. C., Nolasco, R., Rocha, A., and Dubert, J. (2013). Assessing future climate change in the Iberian Upwelling System. *J. Coast. Res.* 65, 1909-1914. doi: 10.2112/S165-323.1

Planque, B., and Arneberg, P. (2017). Principal component analyses for integrated ecosystem assessments may primarily reflect methodological artefacts. *ICES J. Mar. Sci.* 75, 1021-1028. doi: 10.1093/icesjms/fsx223

Polovina, J. J. (2005). Climate variation, regime shifts, and implications for sustainable fisheries. *Bull. Mar. Sci.* 76, 233-244.

R Core Team (2020). *R: A Language and Environment for Statistical Computing*. Vienna: R Foundation for Statistical Computing.

Rochet, M. J., and Trenkel, V. M. (2003). Which community indicators can measure the impact of fishing? A review and proposals. *Can. J. Fish. Aquat. Sci.* 60, 86-99. doi: 10.1139/f02-164

Rodionov, S. N. (2004). A sequential algorithm for testing climate regime shifts. *Geophys. Res. Lett.* 31:L09204. doi: 10.1029/2004GL019448

Rodionov, S. N. (2006). Use of prewhitening in climate change regime shift detection. *Geophys. Res. Lett.* 33:L12707. doi: 10.1029/2006GL025904

Rodionov, S. N., and Overland, J. E. (2005). Application of a sequential regime shift detection method to the Bering Sea ecosystem. *ICES J. Mar. Sci.* 62, 328-332. doi: 10.1016/j.icesjms.2005.01.013

Santos, A. M. P., Borges, M. F., and Groom, S. (2001). Sardine and horse mackerel recruitment and upwelling off Portugal. *ICES J. Mar. Sci.* 58, 589-596. doi: 10.1006/jmsc.2001.1060

Santos, F., Gomez Gesteira, M., and de Castro, M. (2011). Coastal and oceanic SST variability along the western Iberian Peninsula. *Cont. Shelf Res.* 31, 2012-2017. doi: 10.1016/j.csr.2011.10.005

Santos, M. B., González-Quirós, R., Riveiro, I., Cabanas, J. M., Porteiro, C., and Pierce, G. J. (2012). Cycles, trends, and residual variation in the Iberian sardine (*Sardina pilchardus*) recruitment series and their relationship with the environment. *ICES J. Mar. Sci.* 69, 739-750. doi: 10.1093/icesjms/fsr186

Scheffer, M., Carpenter, S. R., Foley, J. A., Folke, C., and Walker, B. H. (2001). Catastrophic shifts in ecosystems. *Nature* 413, 591-596. doi: 10.1038/35098000

Sguotti, C., and Cormon, X. (2018). "Regime shifts - a global challenge for the sustainable use of our marine resources," in *YOU/MARES 8 - Oceans Across Boundaries: Learning from each other*, eds S. Jungblut, V. Liebich, and M. Bode (Cham: Springer), doi: 10.1007/978-3-319-93284-2_11

Shannon, L. J., Christensen, V., and Walters, C. J. (2004). Modelling stock dynamics in the southern Benguela ecosystem for the period 1978-2002. *Afr. J Mar. Sci.* 26, 179-196. doi: 10.2989/18142320409504056

Shannon, L. J., Coll, M., and Neira, S. (2009). Exploring the dynamics of ecological indicators using food web models fitted to time series of abundance and catch data. *Ecol. Indic.* 9:1095. doi: 10.1016/j.ecolind.2008.12.007

Strinimann, L., Conversi, A., and Marini, S. (2019). Detection of regime shifts in the environment: testing "STARS" using systematic and observed time series. *ICES J. Mar. Sci.* 76, 2286-2296. doi: 10.1093/icesjms/fsz148

Sydeman, W. J., Garcia-Reyes, M., Schoeman, D. S., Rykaczewski, R. R., Thompson, S. A., Black, B. A., et al. (2014). Climate change and wind intensification in coastal upwelling ecosystems. *Science* 345, 77-80. doi: 10.1126/science.1251635

Szalaj, D., Torres, M. A., Veiga-Malta, T., Angélico, M. M., Sobrinho-Gonçalves, L., and Chaves, C. (2021). Food-web dynamics in the Portuguese continental shelf ecosystem between 1986 and 2017: unravelling drivers of sardine decline. *Estuar. Coast. Shelf Sci.* 251:107259. doi: 10.1016/j.ecss.2021.107259

Tanner, S. E., Giacomello, E., Menezes, G. M., Mirasole, A., Neves, J., Sequeira, V., et al. (2020). Marine regime shifts impact synchrony of deep-sea fish growth in the northeast Atlantic. *Oikos* 129, 1781-1794. doi: 10.1111/oik.07332
Tomczak, M. T., Heymans, J. J., Yletyinen, J., Niiranen, S., Otto, S. A., and Blenckner, T. (2013). Ecological network indicators of ecosystem status and change in the Baltic Sea. *PLoS One* 8:e75439. doi: 10.1371/journal.pone.0075439

Trigo, R. M., Pozo-Vazquez, D., Osborn, T. J., Castro-Diez, Y., Gamiz-Fortis, S., and Esteban-Parra, M. J. (2004). North Atlantic Oscillation influence on precipitation, river flow and water resources in the Iberian Peninsula. *Int. J. Climatol.* 24, 925–944. doi: 10.1002/joc.1048

Ulanowicz, R. E. (1986). *Growth and Development: Ecosystems Phenomenology*, 2nd Edn. Lincoln, NE: To Excel Press.

Ulanowicz, R. E., and Norden, J. S. (1990). Symmetrical overhead in flow networks. *Int. J. Syst. Sci.* 21, 429–437. doi: 10.1080/00207729008910372

Ullah, H., Leitão, F., Baptista, V., and Chícharo, L. (2012). An analysis of the impacts of climatic variability and hydrology on the coastal fisheries, Engraulis encrasicolus and Sepia officinalis, of Portugal. *Ecohydrol. Hydrobiol.* 12, 337–352. doi: 10.2478/v10104-012-0026-y

Vasconcellos, M., Mackinson, S., Sloman, K., and Pauly, D. (1997). The stability of trophic mass-balance models of marine ecosystems: a comparative analysis. *Ecol. Model.* 100, 125–134. doi: 10.1016/S0304-3800(97)00150-6

Veiga-Malta, T., Szalaj, D., Angélico, M. M., Azevedo, M., Farias, I., Garrido, S., et al. (2019). First representation of the trophic structure and functioning of the Portuguese continental shelf ecosystem: insights into the role of sardine. *Mar. Ecol. Prog. Ser.* 617-618, 323–340. doi: 10.3354/meps12724

Walker, B., Holling, C. S., Carpenter, S. R., and Kinzig, A. (2004). Resilience, adaptability and transformability in social-ecological systems. *Ecol. Soc.* 9:5. doi: 10.5751/ES-00650-090205

Walters, C., Christensen, V., and Pauly, D. (1997). Structuring dynamic models of exploited ecosystems from trophic mass-balance assessments. *Rev. Fish. Biol. Fisher.* 7, 139–172. doi: 10.1023/A:1018479526149

Walters, C. J., Martell, S. J. D., and Korman, J. (2006). A stochastic approach to stock reduction analysis. *Can. J. Fish. Aquat. Sci.* 63, 212–223. doi: 10.1139/f05-213

Weijerman, M., Lindeboom, H., and Zuur, A. F. (2005). Regime shifts in marine ecosystems of the North Sea and Wadden Sea. *Mar. Ecol. Prog. Ser.* 298, 21–39. doi: 10.3354/meps298021

Wooster, W. S., Bakun, A., and McLain, R. M. (1976). Seasonal upwelling cycle along the Eastern boundary of the North Atlantic. *J. Mar. Res.* 34, 131–141.

Wooster, W. S., and Zhang, C. I. (2004). Regime shifts of the North Pacific early indications of the 1976–1977 event. *Prog. Oceanogr.* 60, 183–200. doi: 10.1016/j.pcean.2004.02.005

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

Copyright © 2021 Szalaj, Silva, Ré and Cabral. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.