# Tracking changes in stability of North Sea Atlantic cod over 40 years

Hsiao-Hang Tao, Chih-Hao Hsieh, Manuel Hidalgo, Vasilis Dakos

► To cite this version:

Hsiao-Hang Tao, Chih-Hao Hsieh, Manuel Hidalgo, Vasilis Dakos. Tracking changes in stability of North Sea Atlantic cod over 40 years. ICES Journal of Marine Science, 2025, 82, 10.1093/icesjms/fsaf117. hal-05312389

HAL Id: hal-05312389

https://hal.umontpellier.fr/hal-05312389v1

Submitted on 13 Oct 2025

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
# Tracking changes in stability of North Sea Atlantic cod over 40 years

Hsiao-Hang Tao ID 1,*, Chih-hao Hsieh ID 2,3,4,5, Manuel Hidalgo ID 6, Vasilis Dakos 1

## Abstract

Tracking changes in the stability of low-abundance populations is crucial for effective conservation. While population stability is often evaluated using linear models, many natural populations exhibit state-dependent dynamics, such as the iconic Atlantic cod (Gadus morhua) in the North Sea. This population experienced an abrupt decline around year 2000 and has not shown a clear sign of recovery, despite a reduction in fishing pressure. To investigate the mechanisms underlying its recovery dynamics, we applied a dynamical systems approach to reconstruct the population's state space using yearly abundance data of different age groups between 1982 and 2021. We estimated yearly interaction matrices of age groups and derived a measure of population dynamical stability, defined as the rate at which the population recovers to or departs from its current state following small perturbations to the matrix. Our analysis revealed that the North Sea cod shifted from a disequilibrium state to a stable, low-abundance state in 2003. During this stable phase, some age groups exhibited weaker self-regulation, and age-1 group became more sensitive to fluctuations in the abundance of older groups and showed weakened growth to the next age. These changes likely contributed to the population's failure to recover, even under decreased fishing pressure. The analytical framework used in this study is broadly applicable to assess the stability of structured non-equilibrium populations across different taxa and ecosystems.

**Keywords:** age-specific abundance; non-linear dynamics; empirical dynamical modelling; interaction Jacobians; population stability; recovery dynamics

## Introduction

Many marine populations continue to experience delayed or uncertain rebuilding (Neubauer et al. 2013, Britten et al. 2017, Duarte et al. 2020). Understanding their capacity to withstand disturbance is crucial for sustaining them in challenging environments. The resilience of a fish population is often evaluated by measuring the distance to, or rate of return to, fixed biological reference points, to indicate the potential of recovery (e.g. Neubauer et al. 2013, Trochta et al. 2020). However, some marine fish populations exhibit state-dependent dynamics (i.e. nonlinear dynamics), where population state changes continuously rather than fluctuating around a fixed equilibrium (Anderson et al. 2008, Glaser et al. 2014). For such populations, resilience is more accurately described in relation to their current state at each point in time, rather than to static benchmarks. Tracking how resilience changes over time in these dynamic systems offers valuable insights into the mechanisms underlying recovery dynamics.

Recently, dynamical stability—an indicator based on state-dependent dynamics—was developed to quantify the time-dependent resilience of non-equilibrium natural communities (Ushio et al. 2018). This method uses species time series to reconstruct a multivariate state space based on Takens' embedding theorem (Takens 1981), where the system's trajectory

reflects its intrinsic dynamics. These trajectories are used to estimate time-varying species interactions (Sugihara et al. 2012), and calculate dynamical stability. At any point, the system may be in a stable, unstable, or critical transitional state. In a stable state, it tends to return to its trajectory after perturbations, implying that the state variables (e.g. abundance or biomass) are predictable and likely to stay close to the current condition. In contrast, in an unstable state, the system tends to diverge from its original trajectory after perturbations, signaling the potential for abrupt changes in the state variables. Natural systems often shift between two states, reflecting a continuous state of disequilibrium, where the system alternates between converging toward and diverging from its current trajectory.

While the dynamical stability approach has mainly been applied at the community level, it can also be used to assess population-level stability by reconstructing state space from internal components such as age classes, size groups, or life stages—though this application remains largely unexplored. In marine fish populations, estimating time-dependent interactions between age groups can reveal how these relationships shift over time. Key age group interactions of marine fish populations include growth, where one age group positively influences the next (Hidalgo et al. 2012); reproduc-

<footnote n="1">¹Institut des Sciences de l'Evolution de Montpellier, Université de Montpellier, Montpellier 34095, France</footnote>

<footnote n="2">²Institute of Oceanography, National Taiwan University, Taipei 10617, Taiwan</footnote>

<footnote n="3">³Department of Life Science, National Taiwan University, Taipei 10617, Taiwan</footnote>

<footnote n="4">⁴Institute of Ecology and Evolutionary Biology, National Taiwan University, Taipei 10617, Taiwan</footnote>

<footnote n="5">⁵Research Center for Environmental Changes, Academia Sinica, Taipei 11529, Taiwan</footnote>

<footnote n="6">⁶Spanish Institute of Oceanography (IEO, CSIC), Balearic Oceanographic Center, Group of Ecosystem Oceanography (GRECO), Moll de Ponent s/n, 07190 Palma, Spain</footnote>

<footnote>*Corresponding author. Institut des Sciences de l'Evolution de Montpellier, Université de Montpellier, Montpellier 34095, France E-mail: hsiaohang.tao@gmail.com</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Illustration of the analytical approach in this study. First, the optimal embedding dimension was determined for each age group and for the whole population (Step 1). Then, the nonlinearity of the North Sea cod population was examined (Step 2). After reconstructing the state space (Step 3), we obtained an estimated age group interaction matrix for each year. From each matrix, we performed eigen analysis to extract the absolute value of the dominant eigenvalues as the dynamical stability of that year (Step 4), as well as the absolute value of dominant eigenvector of each age group, standardized by the square root of sum of squares of all eigenvectors of that year to indicate the relative sensitivity to perturbations (Step 5).</figcaption></figure>

tion, with adults contributes to juvenile abundance (Lowerre-Barbieri et al. 2017); cannibalism, where older individuals negatively impact younger ones (Pereira et al. 2017); and self-regulation, where competition within age groups generates negative feedbacks (Moksnes 2004). These interactions are often abundance-dependent: reduced adult numbers may weaken reproduction and juvenile recruitment, while reduced overall abundance can ease self-regulation within ages groups and enhance survival of the remaining individuals.

Using state-space reconstruction of age groups, we can potentially assess overall population stability and also identify which age groups are more sensitive to fluctuations in abundance, by comparing the relative change in each age group's distance from its original trajectory following a disturbance—the greater the deviation, the higher the sensitivity of an age group. This method has been developed at the community level, to rank relative sensitivity of species to perturbations of species (Medeiros et al. 2022), and it is potentially applicable to the age-structured populations.

Atlantic cod in the North Sea is an iconic example of a fish population experiencing delayed recovery. As a target commodity in one of the world's oldest fisheries, North Sea cod have shaped the history, economy, and culture of Western Europe since the 16th century (Goodlad 2022). However, this population now faces multiple challenges, including rising temperatures, more frequent and intense extreme events (IPCC 2023), and the risk of irreversibility regime shifts in the North Sea ecosystem (Sguotti et al. 2022). Over recent decades, the population has undergone several abrupt declines, with the most significant around 2000 (Blöcker et al. 2022), primarily due to a substantial loss of juveniles (ICES 2024). Year 2000 also marked the start of reduced fishing mortality, as a consequence of reforms of the European Union's Common Fisheries Policy (European Parliament 2013, Salomon et al. 2014). Despite eased fishing pressure, the population size remains low. This raises to a critical question: whether and how do changes in dynamical stability and interactions among age group over time may explain of the population dynamics of North Sea cod?

Using survey data of North Sea cod from 1982 to 2021, this study addresses the following question: How have dy-

namical stability, interactions among age groups, and the relative sensitivity of each age group to perturbations changed over time? We hypothesize that North Sea cod entered a low-abundance stable state after 2000 due to altered age group interactions and imbalanced sensitivity among age groups, causing the population to remain at low-abundance level despite relaxed fishing pressure.

## Materials and methods

### Data

We obtained the abundance data of six age group of the North Sea cod between 1982 and 2021 from the International Bottom Trawl Survey (IBTS) of the International Council for the Exploration of the Sea (ICES) (https://data.ices.dk). We used the yearly winter data (from January to March) because it has higher sampling consistency than other quarters. Abundance was reported as catch per unit effort (CPUE), measured as the number of fish caught per hour of trawling within each survey rectangle (1° longitude by 0.5° latitude) across the North Sea. To calculate the abundance for each age group in the entire North Sea, we summed the CPUE data from all rectangles. To detect potential abrupt changes in total population abundance over time, we applied a recently developed change point analysis (Pélissié et al. 2024). We obtained yearly estimates of mean fishing mortality for age 2–4 groups from the ICES stock assessment report (ICES 2020).

### Analytical framework

Using abundance time series of six age groups, we followed a five-step analytical workflow (Fig. 1). First, we determined the optimal embedding dimension for each age group and for the whole population (Step 1). Next, we assessed the nonlinearity of the North Sea cod population dynamics to justify the use of dynamical systems method (Step 2). We then estimated the time-dependent, pairwise age group interactions to obtain interaction matrices (Step 3). From these matrices, we performed eigen analysis to extract the absolute value of the dominant eigenvalues as a proxy for population's dynamical stability (Step 4), as well as the absolute value of the dominant eigenvector (after standardization) for each age group to as-
sess the relative sensitivity of each age group to perturbations (Step 5).

## Step 1: Determine the optimal embedding dimension

We first used simplex projection (Sugihara and May 1990) to determine the optimal embedding dimension for the whole population and for each age group. The embedding dimension refers to the optimal number of state variables, or the optimal number of time lags of a state variable. This information is needed to reconstruct the manifold in the state space to test the nonlinearity of the population (Step 2) and construct interaction matrix of age groups (Step 3).

To begin with, we divided the time series of the population or an age group into a library set and a prediction set. Within the library set, we created state-space reconstructions for different embedding dimensions ranging from 2 to 10. For each dimension, we constructed geometric polygons, or simplex, from neighboring points around a given point. These simplex were then used to predict the future states of the system. Specifically, for each point in the state space, the future state was predicted based on the weighted average of the future states of the points that form the simplex. The weights were inversely related to the distances from the points to the vertex of the simplex.

We then determined the optimal embedding dimension by comparing and choosing the dimension that provides the best prediction skill. The highest prediction skill is referred by the lowest root mean square error (RMSE) between the predicted states to the observed states in the prediction set. The resulting optimal embedding dimension for the whole population is 5 and for age groups 1 to 6 are 4, 2, 8, 10, 10, 10, respectively.

## Step 2: Identify the nonlinearity of the population

We used the S-map method (Sugihara 1994, Hsieh et al. 2005) to assess whether the North Sea cod population exhibits nonlinear dynamics. Based on Step 1, the optimal embedding dimension for the overall population was determined to be 5. Accordingly, we reconstructed the system's state space using the population abundance time series and four time-lagged variables. For each point in this reconstructed state space, a locally weighted linear regression is applied using all available points to predict future states of the system. The degree of nonlinearity is governed by a tuning parameter, $\theta$, which controls the weighting of neighboring points. When $\theta = 0$, all neighbors are weighted equally, resulting in a purely linear model. As $\theta$ increases above 0, closer neighbors receive greater weight, allowing the model to capture more localized and nonlinear dynamics. Thus, a higher $\theta$ indicates stronger nonlinearity in the population dynamics.

To determine the optimal value of $\theta$, we ran the S-map analysis using a range of $\theta$ values from 0 to 8 and selected the one that produced the best forecast skill, as indicated by the lowest mean absolute error (MAE). Then, we followed a procedure described by Clark and Luis (2019) to evaluate the significance of nonlinearity for each age group. We first calculated ΔMAE as the difference between MAE at $\theta = 0$ and the minimum MAE. Next, to generate null distribution, we used phase-randomized surrogates to create 1000 abundance time series. These time series preserve the basic statistical properties of the time series such as autocorrelation. We ran the S-map analysis with each of these time series and computed ΔMAE from the null distribution. The dynamics are considered sig-

nificantly nonlinear if ΔMAE from the original time series is higher than the 97.5% quantile of the null distribution.

## Step 3: Estimate age group interactions

We quantified the strength and the direction of interaction between each pair of age groups through time by reconstructing the state space of the system, using the Multiview Distance Regularized S-map method (MDR S-map) (Chang et al. 2021). This approach reconstructs multiple low-dimensional state spaces with the optimal embedding dimension to determine neighborhood relationships among data points. These distances are then ensembled to determine the weights used in the multivariate S-map.

In this study, we used six age groups to reconstruct the low-dimensional state space, with maximum embedding dimension of six. For age groups 3–6, where the optimal embedding dimensions exceed six, we limited them to six. The interaction strength of age group *i* at year *t* on age group *j* in the following year is approximated as the interaction Jacobian ($\partial A_{i(t+1)}/\partial A_{i(t)}$), where $A_{i(t)}$ is the abundance of age *i* at year *t*. For self-interactions, defined as the change of any age group influences on itself ($\partial A_{i(t+1)}/\partial A_{i(t)}$), we subtract the value of interaction Jacobian by 1 (Miki et al. 2024). We compiled all the pairwise interaction Jacobians between the six age groups to generate a 6 × 6 discrete-time Jacobian matrix that is interpreted as the interaction network of our six age groups.

The interaction matrices are estimated based on the distance between time points *t* and *t* + 1 in the state space. As a result, the time series of age group interactions is one year shorter, with 39 time points from 1982 to 2020. This adjustment also applies to the two indices derived from the interaction matrix: the population stability and the relative sensitivity of age groups.

## Step 4: Compute the population stability

For each year, we applied eigen analysis to the age group interaction matrix to calculate the absolute value of the dominant eigenvalue (|DEV|), following Ushio et al. (2018), as a proxy for the population's dynamical stability. This metric reflects "whether system trajectories at a given time are converging or diverging" (Ushio et al. 2018). Specifically, it quantifies how the system responds to small perturbations—whether it tends to return to its original trajectory or deviate further—without requiring the system to be at equilibrium. Thus, stability is assessed locally, based on the system's state at each time point, rather than relative to a fixed equilibrium. A |DEV| of 1 indicates the system is at a critical transition point. Values of |DEV| <1 suggest local stability, with smaller values indicating greater stability. Conversely, |DEV| values greater than 1 indicate local instability, with larger values representing increased instability.

## Step 5: Quantify the relative sensitivity of age groups to perturbations

To quantify the relative sensitivity of each age group to perturbations, we used the eigenvector ranking method described by Medeiros et al. (2022). As in Step 4, we performed eigen analysis on the annual interaction matrix of age groups. We then extracted the dominant eigenvector for each age group and each year, and standardized the value by dividing each by the square root of the sum of squares of all six eigenvectors for that year. The absolute value of each standardized eigenvector represents the relative sensitivity of the corresponding
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Changes in the abundance and stability of the North Sea cod population over time. (a) Winter abundance of the population. Vertical lines indicate break points detected in this study and previous work: 1996 (Blöcker et al. 2023), 1999 (Blöcker et al. 2023), 2000 (this study), 2006 (Blöcker et al. 2022), and 2007 (Blöcker et al. 2023). Year 2000 also marks a beginning of decline in fishing mortality after the fisheries policy reform (Salomon et al. 2014). (b) Abundance of each age group as catch per unit effort (CPUE). (c) Population dynamical stability as the absolute value of the dominant eigenvalue (|DEV|) derived from reconstructed age group interaction matrices. Dotted horizontal line indicates the critical transition at |DEV| = 1. A population is considered stable when |DEV| < 1 and unstable when |DEV| > 1. The color indicates three phases based on the fluctuation pattern of |DEV|. (d) Mean fishing mortality between age 2 and 4 of the population.</figcaption></figure>

age group to perturbations. We also ranked these values from 1 to 6 to compare sensitivity levels across age groups. This indicator approximates the distance of each age group deviates from its unperturbed trajectory following a small perturbation or fluctuations in abundances. Age groups with higher sensitivity values are more likely to diverge from their original trajectory, indicating that their abundance is strongly influenced by changes in the abundances of other age groups.

## Results

### Population dynamics and fishing pressure of North Sea cod

Between 1982 and 2021, we identified an abrupt decline in the total abundance of North Sea cod in the year 2000 using change point analysis (Fig. 2a, Supplementary Fig. S1).

Previous studies based on biomass time series also reported regime shifts in 1996, 1999, 2006, and 2007 (Blöcker et al. 2022, Blöcker et al. 2023) (Fig. 2a). Among the six age groups, age-1 and age-2 cod consistently dominated the population throughout the study period. (Fig. 2b). After 2000, the total abundance of the population dropped by >60%, largely driven by sharp declines in these two younger age groups. Around the same time, the mean fishing mortality for age-2 to age-4 cod began to decrease (Fig. 2d) as a consequence of fishing policy reforms (European Parliament 2013, Salomon et al. 2014).

### Stability of North Sea cod

The North Sea cod population exhibits nonlinear dynamics (Supplementary Fig. S2), supporting the use of the dynamic
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Age group interaction network of 1982–2000 and 2001–2021. The circles indicate age groups, and the arrows indicate the mean interaction strength for each pair of age groups averaged from yearly values within period of 1982–2000 (a) and 2001–2020 (b). The thickness of the arrows indicates the interaction strength. To facilitate visualization, the mean interaction strength was classified into four categories: weak, medium, strong, and very strong.</figcaption></figure>

stability index (|DEV|) to assess its resilience to perturbations. We identified three distinct phases of the study period based on the temporal fluctuations of |DEV|: (1) From 1982 to 2002, |DEV| fluctuated between values above and below 1, indicating alternating periods of local stability and instability, and indicating that the population was at disequilibrium state. (2) Between 2003 and 2014, |DEV| remained consistently below 1, indicating that the population was at a stable state. (3) From 2015 to 2020, |DEV| again showed fluctuations above and below 1, indicating a return to a disequilibrium state (Fig. 2c).

## Age group interactions

We compared the mean interaction strength of each age group pair between two periods: 1982–2000 and 2001–2020 (Fig. 3, Supplementary Fig. S3). The facilitate interpretation, we categorized the interaction strength based on the absolute value of jacobian interaction: weak (<0.3), medium (0.3–0.6), strong (0.6–0.9), and very strong (>0.9).

In both periods, the strongest interactions were self-interactions within age groups—that is, how changes in an age group's abundance affect its own abundance at the next time point. These self-interactions were consistently negative, with mean strengths ranging from medium to very strong. The second most prominent interactions were growth, describing how changes in one age group influence the next age group at the following time step. These were consistently positive, with

mean strengths from medium to strong. Additionally, a strong positive link from the age-2 to age-1 group was observed in 2001–2020. All other interactions in both periods were classified as weak.

The most notable change between the two periods was a weakening of the negative self-interactions in the age-1, age-2, and age-4 groups, which shifted from very strong to strong. In contrast, the negative self-interaction within the age-6+ group strengthened, increasing from medium to strong. Additionally, the positive growth interaction from age-1 to age-2 weakened from strong to medium, while the positive interaction from age-2 to age-1 strengthened markedly from weak to strong.

The results of age group interactions were mostly identical when comparing the interactions between two first phases of dynamical stability (1982–2002 and 2003–2014) (Supplementary Fig. S4), except that the positive link from age-2 to age-1 strengthened from weak to medium, and that the presence or sign of few weak interactions differed.

## Relative sensitivity of each age group to perturbations

We examined whether and how the relative sensitivity of six age groups has changed over the study period by comparing the absolute value of standardized dominant eigenvector between each age group. We used violin plots to visualize the
yearly variability in the relative sensitivity for each age group within periods of 1982–2000 and 2001–2020 (Fig. 4).

Comparing two periods, the relative sensitivity of age
groups exhibited higher yearly variability in the period of
1982–2000. This indicates greater year-to-year fluctuations
in the relative sensitivity of age groups. In contrast, the rel-
ative sensitivity of age groups exhibited lower yearly variabil-
ity in the period of 2001–2020. This suggests more consis-
tent sensitivity values across years. Particularly, age-1 group
exhibited the highest sensitivity for most years. This find-
ing was robust when using the ranking of this index, al-
though with greater year-to-year variation for two peri-
ods (Supplementary Fig. S5). The results were also identi-
cal when comparing between two first phases of dynam-
ical stability (1982–2002 and 2003–2014) (Supplementary
Fig. S6).

Discussion

Although fishing pressure on many fisheries has decreased in recent decades, more than half of the assessed marine fish populations worldwide still require rebuilding (Worm et al. 2009). Understanding how resilient these populations are to disturbances can help explain their population dynamics. In this study, we focused on the North Sea cod population, which experienced a sharp decline in 2000. Using a dynamical systems approach, we analyzed the population's stability, the interactions among different age groups, and the sensitivity of each age group to perturbations, and examined how these properties have changed over time.

Population stability and dynamics of North Sea cod

We identified three distinct phases in the dynamical stability of the North Sea cod over the study period. During the first phase between 1982 and 2002, the population was at a dis-equilibrium state, implying high unpredictability and a ten-dency for abrupt shift in the population state. This pattern is consistent with the observed fluctuations in abundance dur-ing this phase, followed by a sharp decline around the year 2000. During the second phase between 2003–2014, the pop-ulation entered a stable state, indicating greater predictability and consistent abundance levels. This finding corresponds to consistently low abundance throughout this phase despite of reduced fishing pressure from 2000, implying a lack of capac-ity to recovery. During the third phase between 2015–2020, the population returned to a disequilibrium state, reflected by the fluctuating abundance during this phase, with peaks in 2015 and 2017 followed by a decline in 2019. By 2020, the population remained in a disequilibrium state, suggesting a potential change in the state but the direction is unpredictable. In the best scenario, the population could transition to a stable state with higher abundance, signaling recovery. Alternatively, it may shift to a stable state with even lower abundance, indi-cating a further decline in population size.

It is intuitive to associate a stable population state with a
healthy condition and an unstable state with an unhealthy
one. However, by definition, dynamical stability does not di-
rectly reflect the health of a population. Rather, it describes
the system's tendency to either remain near its current state
or shift away from it—that is, whether it returns to its
original trajectory after a disturbance or continues to devi-
ate from it (Ushio et al. 2018). Assessing population health,

therefore, requires considering both the population dynamics
and its stability. For example, a high-abundance stable state
may indicate a healthy population, as it tends to returns to
high-abundance level after perturbations. In contrast, a low-
abundance stable state may signal risk, as the population is
resistant to change and may be "trapped" in a depleted con-
dition. This was the case for the North Sea cod during 2003–
2014. Alternatively, when a population is at an unstable state,
or frequent transitions between stable and unstable states, in-
dicates that it is undergoing change, such as increases or de-
creases in abundance or biomass. During 1982–2002, the dis-
equilibrium state of the North Sea cod likely contributed to
the abrupt decline observed around 2000. From 2015 to 2020,
the population was again unstable, and it remains uncertain
whether it will transition toward recovery or further decline.

Based on our findings and previous literature, the North
Sea cod likely experienced one or multiple regime shifts
between 1996 and 2007. The exact timing of these shifts
varies depending on the method used and state variable an-
alyzed (e.g. biomass or abundance) used. Chronologically,
Blöcker et al. (2023) detected a shift in 1996 based on a non-
stationary functional relationship between spawning stock
biomass (SSB) and recruitment, followed by a shift in 1999
identified through an abrupt change in the SSB. Our study sug-
gests a shift in 2000 based on an abrupt change in abundance,
and another in 2003 inferred from changes in dynamical sta-
bility patterns. Further, Blöcker et al. (2022) reported a shift in
2006 based on an abrupt change in SSB, and finally, Blöcker et
al. (2023) identified a regime shift in 2007 through evidence
of hysteresis in the relationship between SSB and fishing mor-
tality.

Age group interactions and relative sensitivity of age groups to perturbations

To further understand the mechanisms underlying population
dynamics, we analyzed changes in age group interactions and
the relative sensitivity of each age group to perturbations over
time. Specifically, we compared two timeframes: before and
after the abrupt decline in abundance (1982–2000 vs. 2001–
2020), and the two initial phases of dynamical stability (1982–
2002 vs. 2003–2014). Our analysis revealed that, in the latter
period, the self-interaction strengths of juveniles (age-1 and
age-2 groups) decreased. This pattern may reflect a positive
density-dependent effect (Hilborn et al. 2014, Perälä et al.
2022). The sharp decline in juvenile abundance may have re-
duced cannibalism and intra-group competition, thereby re-
laxing self-regulation and enhancing juvenile survival.

However, despite the relaxed self-regulation observed in
younger age groups, the strength of growth from one age
to the next did not increase in the latter period. In fact,
the strength of growth from age-1 to age-2 group declined.
This suggests that juvenile cod survival is heavily influenced
by external pressures such as warming (Perry et al. 2024),
fishing—particularly targeting age 2–4 groups (Möllmann et
al. 2021)—and predations (Swain and Sinclair 2000). These
effects are particularly pronounced when the population de-
mography is eroded (Anderson et al. 2008, Hsieh et al. 2010,
Hidalgo et al. 2012). Additionally, the age-1 group has be-
come more sensitive to fluctuations in the abundance of older
groups, likely reflecting impaired reproduction due to the
extremely low number of mature adults. In summary, these
findings suggest that while positive density-dependent pro-
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Distribution of the yearly relative sensitivity to perturbation of each age group for 1982–2000 and 2001–2020. Each point represents the relative sensitivity of an age group compared to other age groups in a given year, computed as the absolute value of dominant eigenvector derived from the age group interaction matrix. Each value was standardized by the square root of the sum of squares of all six eigenvectors for a given year. The distribution of yearly relative sensitivity for each age group is summarized using boxplots and violin plots for 1982–2000 (a) and the 2001–2020 (b).</figcaption></figure>

cesses may have supported population growth, the high vul-
nerability to external pressures of age-1 group and low repro-
ductive output from declined number of adults likely limited
the survival of the age-1 group. These processes may have in-
creased the sensitivity of age-1 group to fluctuations in abun-
dance of other age groups, and weakened its growth into sub-
sequent age group, which ultimately dampened the entire life
cycle. These mechanisms provide a potential explanation on
why North Sea cod have remained in a low-abundance state
since 2003, despite reductions in fishing pressure.

Limitations

This study has several limitations. First, we assessed the stabil-
ity of the North Sea cod population as a whole, even though
recent evidence suggests it functions as a metapopulation com-
posed of multiple subpopulations (Holmes et al. 2014). Future
research could explore differences in stability and age group
interactions among these subpopulations to gain a deeper un-
derstanding of metapopulation dynamics. Second, it is impor-
tant to note that our findings are based on winter data, while
interactions or population stability may vary seasonally due to
behavioral changes. For instance, during winter, mature cod
undertake spawning migrations and become more spatially
aggregated than in other seasons (Neat et al. 2014), which
may influence interactions among adults and between adults
and juveniles. Third, our analysis focused on strong age group
interactions—such as growth and self-regulation—rather than
weak interactions, which also contribute to population re-
silience. However, many of the weak links identified in this

study were difficult to interpret and may reflect random noise.
At present, there is a lack of method for evaluating the signif-
icance of such interactions, highlighting the need for further
methodological development.

Methodology consideration

There are several considerations regarding the use of the dynamical stability approach. First, dynamical stability values tend to be underestimated due to process and measurement errors (Grziwotz et al. 2023). There is currently no available method for quantifying the magnitude of these errors. Second, reliable reconstruction of the state space typically requires at least 35 data points (Sugihara et al. 2012), with accuracy improving when longer and less noisy time series are available (Chang et al. 2017). These constraints reveal the uncertainty of the findings and interpretations. Finally, the dynamical stability index does not provide information on the direction of change, nor does it allow for quantification of the probability of a shift. For predictive purposes, alternative approaches such as early warning signals may be more suitable (e.g. Clements et al. 2019, Grziwotz et al. 2023).

There are important differences between the dynamical sta-
bility method used in this study and the approaches commonly
found in the regime shift literature. Studies on regime shifts
(e.g. Sguotti et al. 2019, Möllmann et al. 2021, Blöcker et al.
2023) typically identify the timing of regime changes by ana-
lyzing how the state of a population or community responds
to external drivers. For example, detecting hysteresis or bi-
furcations in biomass in response to factors such as fishing
mortality or temperature. In contrast, the dynamical stability method focuses on whether system trajectories are converging or diverging at each point in time. It does not identify or quantify basins of attraction, nor does it require the system to be at equilibrium. This method detects whether perturbations are amplified or dampened by the interacting components within the system (e.g. age or size group interactions within a population). Although the two approaches differ fundamentally in their focus and assumptions, both offer valuable insights. Integrating findings across diverse methods can enhance our understanding of population status and dynamics.

## Perspectives and management implications

The lack of recovery in the North Sea cod population observed in this study aligns with broader patterns seen in marine fish populations, where prolonged recovery is common following long histories of exploitation (Neubauer et al. 2013). Our findings highlight the critical role of improving survival in the age-1 group, which is both highly sensitive to perturbations and shows impaired transition to the next age. One potentially effective measure is reducing bycatch and discard rates of age-1 groups. This is particularly relevant given that commercial catches in this fishery are dominated by juveniles under three years old, with bycatch and discards accounting for an estimated 23% of the total catch by weight in 2023 (ICES 2024). In addition, further reductions in fishing mortality on adult cod are essential to prevent continued declines in adult abundance and to support reproductive output and recruitment (Blöcker et al. 2022). Together, these measures could hopefully contribute to rebuilding the population and enhancing its resilience to environmental pressures.

## Acknowledgments

We thank the research and survey teams of International Bottom Trawl Survey Working Group from International Council for the Exploration of the Sea (ICES) for providing the data. We thank Alejandro Cano, Benoît Pichon, and Ismaël Lajaaiti for commenting on an earlier version. We thank Masayuki Ushio for insightful discussion.

## Author contributions

Hsiao-Hang Tao (Conceptualization [lead], Formal analysis [lead], Methodology [lead], Writing – original draft [lead]), Chih-hao Hsieh (Conceptualization [equal], Methodology [equal], Validation [equal], Writing – review & editing [equal]), Manuel Hidalgo (Conceptualization [equal], Validation [equal], Writing – review & editing [equal]), and Vasilis Dakos (Conceptualization [equal], Methodology [equal], Project administration [equal], Supervision [equal], Validation [equal], Writing – review & editing [equal])

## Supplementary data

Supplementary material is available at the *ICES Journal of Marine Science* online version of the manuscript.

*Conflict of interest:* All authors declare that they have no conflicts of interest.

## Funding

This work is supported by Marie Curie postdoctoral fellowship of H.H.-T. (HORIZON-MSCA-2021-PF-01, European Commission, LSP 244267; MAFIS project). C.H.-H. is supported by the National Science and Technology Council, Taiwan. M.H. acknowledges CORES project (PID2023-151282OB-I00, Spanish Ministry of Science and Innovation, MCIN/AEI/10.13039/501100011033).

## Data availability

All the data used this study were downloaded from open-access online databases as mentioned above. We used “rEDM” package (https://cran.r-project.org/src/contrib/Archi ve/rEDM, version 1.2.3) to perform all the analysis. The MDR S-map was performed using functions provided in Chang et al. (2021). All statistical analyses were performed in R 4.3.1. Data and codes for reproducing the results are deposited in https://github.com/HsiaoHang/git_north_sea_code_stability.

## References

* Anderson CNK, Hsieh C, Sandin SA *et al.* Why fishing magnifies fluctuations in fish abundance. *Nature* 2008;452:835–9. https://doi.org/10.1038/nature06851

* Blöcker A, Sguotti C, Möllmann C. Discontinuous dynamics in North Sea cod Gadus morhua caused by ecosystem change. *Mar Ecol Prog Ser* 2022;713:133–49. https://doi.org/10.3354/meps14342

* Blöcker AM, Gutte HM, Bender RL *et al.* Regime shift dynamics, tipping points and the success of fisheries management. *Sci Rep* 2023;13:289. https://doi.org/10.1038/s41598-022-27104-y

* Britten GL, Dowd M, Kanary L *et al.* Extended fisheries recovery timelines in a changing environment. *Nat Commun* 2017;8:15325. https://doi.org/10.1038/ncomms15325

* Chang C, Miki T, Ushio M *et al.* Reconstructing large interaction networks from empirical time series data. *Ecol Lett* 2021;24:2763–74. https://doi.org/10.1111/ele.13897

* Chang C-W, Ushio M, Hsieh C. Empirical dynamic modeling for beginners. *Ecol Res* 2017;32:785–96. https://doi.org/10.1007/s11284-017-1469-9

* Clark TJ, Luis AD. Nonlinear population dynamics are ubiquitous in animals. *Nat Ecol Evol* 2019;4:75–81. https://doi.org/10.1038/s4159 5-019-1052-6

* Clements CF, McCarthy MA, Blanchard JL. Early warning signals of recovery in complex systems. *Nat Commun* 2019;10:1681. https://doi.org/10.1038/s41467-019-09684-y

* Duarte CM, Agusti S, Barbier E *et al.* Rebuilding marine life. *Nature* 2020;580:39–51. https://doi.org/10.1038/s41586-020-2146-7

* European Parliament. Regulation (EU) No 1380/2013 of the European Parliament and of the Council. *Off J Eur Union* 2013;354:22–61

* Glaser SM, Fogarty MJ, Liu H *et al.* Complex dynamics may limit prediction in marine fisheries. *Fish Fish* 2014;15:616–33. https://doi.org/10.1011/j.1600-0587.2011.07314.x

* Goodlad J. *The Salt Roads: How Fish Made a Culture*. London: Birlinn, 2022

* Grziwotz F, Chang C-W, Dakos V *et al.* Anticipating the occurrence and type of critical transitions. *Sci Adv* 2023;9:abq4558. https://doi.org/10.1126/sciadv.abq4558

* Hidalgo M, Rouyer T, Bartolino V *et al.* Context-dependent interplays between truncated demographies and climate variation shape the population growth rate of a harvested species. *Ecography* 2012;35:637–49. https://doi.org/10.1111/j.1600-0587.2011.07314.x

* Hilborn R, Hively DJ, Jensen OP *et al.* The dynamics of fish populations at low abundance and prospects for rebuilding and recovery. *ICES J Mar Sci* 2014;71:2141–51. https://doi.org/10.1093/icesjms/fsu035
Holmes SJ, Millar CP, Fryer RJ *et al.* Gadoid dynamics: differing perceptions when contrasting stock vs. population trends and its implications to management. *ICES J Mar Sci* 2014;71:1433-42. https://doi.org/10.1093/icesjms/fsu075

Hsieh C, Glaser SM, Lucas AJ *et al.* Distinguishing random environmental fluctuations from ecological catastrophes for the North Pacific Ocean. *Nature* 2005;435:336-40. https://doi.org/10.1038/natu re03553

Hsieh C, Yamauchi A, Nakazawa T *et al.* Fishing effects on age and spatial structures undermine population stability of fishes. *Aquat Sci* 2010;72:165-78. https://doi.org/10.1007/s00027-009-0 122-2

ICES,ICES Working Group on the Assessments of demersal stocks in the North Sea and Skagerrak (WGNSSK). *ICES Scientific Reports*, 2020. https://doi.org/10.17895/ices.pub.6092 (01 April 2023, date last accessed).

ICES,Working Group on the Assessment of Demersal Stocks in the North Sea and Skagerrak (WGNSSK). *ICES Scientific Reports*,2024. 198122080 Bytes, 1249. https://doi.org/10.17895/ices.pub.256056 39(26 July 2024, date last accessed)

IPCC. *Climate Change 2022-Impacts, Adaptation and Vulnerability: Working Group II Contribution to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change*. 1st ed. Cambridge University Press, 2023.

Lowerre-Barbieri S, DeCelles G, Pepin P *et al.* Reproductive resilience: a paradigm shift in understanding spawner-recruit systems in exploited marine fish. *Fish Fish* 2017;18:285-312. https://doi.org/10.1111/faf.12180

Medeiros LP, Allesina S, Dakos V *et al.* Ranking species based on sensitivity to perturbations under non-equilibrium community dynamics. *Ecol Lett* 2022;26:170-83. https://doi.org/10.1111/ele.14131

Miki T, Chang C-W, Ke P-J *et al.* How to quantify interaction strengths? A critical rethinking of the interaction jacobian and evaluation methods for non-parametric inference in time series analysis. *Physica D: Nonlinear Phenomena* 2024;476:134613. https://www.scienc edirect.com/science/article/pii/S0167278925000922

Moksnes P-O. Self-regulating mechanisms in cannibalistic populations of juvenile shore crabs *Carcinus maenas*. *Ecology* 2004;85:1343-54. https://doi.org/10.1890/02-0750

Möllmann C, Cormon X, Funk S *et al.* Tipping point realized in cod fishery. *Sci Rep* 2021;11:14259. https://doi.org/10.1038/s41598-0 21-93843-z

Neat FC, Bendall V, Berx B *et al.* Movement of Atlantic cod around the British Isles: implications for finer scale stock management. *J Appl Ecol* 2014;51:1564-74. https://doi.org/10.1111/1365-2664.12343

Neubauer P, Jensen OP, Hutchings JA *et al.* Resilience and recovery of overexploited marine populations. *Science* 2013;340:347-9. https://doi.org/10.1126/science.1230441

Pélissié M, Devictor V, Dakos V. A systematic approach for detecting abrupt shifts in ecological timeseries. *Biol Conserv* 2024;290:110429. https://doi.org/10.1016/j.biocon.2023.110429

Perälä T, Hutchings JA, Kuparinen A. Allee effects and the Allee-effect zone in northwest Atlantic cod. *Biol Lett* 2022;18:20210439. https://doi.org/10.1098/rsbl.2021.0439

Pereira LS, Agostinho AA, Winemiller KO. Revisiting cannibalism in fishes. *Rev Fish Biol Fish* 2017;27:499-513. https://doi.org/10.1007/s11160-017-9469-y

Perry D, Tamarit E, Sundell E *et al.* Physiological responses of Atlantic cod to climate change indicate that coastal ecotypes may be better adapted to tolerate ocean stressors. *Sci Rep* 2024;14:12896. https://doi.org/10.1038/s41598-024-62700-0

Salomon M, Markus T, Dross M. Masterstroke or paper tiger—The reform of the EU's Common Fisheries Policy. *Mar Policy* 2014;47:76-84. https://doi.org/10.1016/j.marpol.2014.02.001

Sguotti C, Blöcker AM, Färber L *et al.* Irreversibility of regime shifts in the North Sea. *Front Mar Sci* 2022;9:945204. https://doi.org/10.3891/fmars.2022.945204

Sguotti C, Otto SA, Frelat R *et al.* Catastrophic dynamics limit Atlantic cod recovery. *Proc R Soc B Biol Sci* 2019;286:20182877. https://doi.org/10.1098/rspb.2018.2877

Sugihara G. Nonlinear forecasting for the classification of natural Time Series. *Philos Trans Phys Sci Eng* 1994;348:477-95.

Sugihara G, May R, Ye H *et al.* Detecting causality in complex ecosystems. *Science* 2012;338:496-500. https://doi.org/10.1126/science.1227079

Sugihara G, May RM. Nonlinear forecasting as a way of distinguishing chaos from measurement error in time series. *Nature* 1990;344:734-41. https://doi.org/10.1038/344734a0

Swain DP, Sinclair AF. Pelagic fishes and the cod recruitment dilemma in the Northwest Atlantic. *Can J Fish Aquat Sci* 2000;57:1321-5. https://doi.org/10.1139/f00-104

Takens F. DA Rand, L-S Young. *Dynamical Systems and Turbulence*. New York, Springer 1981:366-81.

Trochta JT, Branch TA, Shelton AO *et al.* The highs and lows of herring: a meta-analysis of patterns and factors in herring collapse and recovery. *Fish Fish* 2020;21:639-62. https://doi.org/10.1111/faf.12452

Ushio M, Hsieh C, Masuda R *et al.* Fluctuating interaction network and time-varying stability of a natural fish community. *Nature* 2018;554:360-3. https://doi.org/10.1038/nature25504

Worm B, Hilborn R, Baum JK *et al.* Rebuilding global fisheries. *Science* 2009;325:578-85. https://doi.org/10.1126/science.1173146

Handling editor: David M Kaplan