We used a previously published time series data set (Eikeset et al., 2016) on the body length of the northeast Arctic cod (Gadus morhua), to study the dynamics of a phenotypic trait. The data set consists of mean body length for different age stages (age 3–12 years) of the population from year 1946 to 2004 (Supplementary Figure 1). The original data set includes body length for all 10 age stages. For simplicity, the data set was divided into three broader age classes: age class I (age 3–5), age class II (age 6–9), and age class III (age 10–12). Each class has a different competition (compete for the shared resources) ability and reproductive (egg supply) rate per biomass. We assume that cod of the age class I has the lowest competition ability (characterized by the maximal attack rate) and reproductive rate, the age class II has a moderate level, and the age class III has a high level; we assume that three age classes have low, moderate, and high growth potential in body length (see Table 1), respectively. In each age class, means of body length were calculated. Finally, to facilitate data analysis, three time series were transformed by a square-root power transformation to suppress sharp peaks.

For the empirical data, we employed a state-space reconstruction (i.e., we reconstructed the multidimensional dynamics of each body length class, using the body length time series) (Takens, 1981; Becks et al., 2005; Benincà et al., 2008; Wang et al., 2019; Dakos, 2020). To this end, the C–C method (Kim et al., 1999), which is useful for smaller data sets, was first used to calculate the time delay (τ) of each time series. Second, combined with the time delay (τ), the Grassberger–Procaccia (G–P) method was used to calculate the embedding dimension (m) (Grassberger and Procaccia, 1983). Third, by combining τ and m, the largest LE of each time series was derived (Wolf et al., 1985); moreover, in order to distinguish whether chaos is driven by external environmental noise or endogenous factors, we use wdencmp function (Donoho et al., 1995) to filter out the environmental noise of empirical data and re-calculate the largest LE of body length (see Supplementary Figure 2 and simulation codes in Supplementary material). Notably, the system has chaotic dynamics only if the Lyapunov exponent is larger than zero. The Lyapunov exponent quantifies the rate of exponential divergence (or convergence) of nearby trajectories (Strogatz, 1994), and a positive LE indicates chaos where the magnitude of LE effectively measures the system sensitivity to initial conditions.

Different with the attractor reconstruction approach, we combine ecological processes (consumer-resource dynamics) with evolutionary processes (phenotypic trait dynamics) to simulate the complex dynamics in fish populations.

The time series data show that fish in the different age stages have different body lengths (Supplementary Figure 1). Different with early age-structured consumer-resource models (Schreiber and Rudolf, 2008; Nilsson et al., 2018), in this study, 10 ages are divided into three age classes (mean, body, and length), under the biologically plausible assumptions (i.e., different maximal attack rate and reproductive rate in Table 1) of each age class, the system is similar to a multi-species consumer-resource system. Here, we establish a bioenergetic consumer-resource model (McCann, 1998) to include phenotypic representations of quantitative trait evolution (Cortez, 2016, 2018; McPeek, 2017; Yamamichi and Letten, 2021). The model can be written as:

Where i = 1, 2, 3, and

The state variable R is the resource biomass; C_i, J_i, and S_i are the biomass, egg biomass and body length of fish in age class i, respectively; S_R is the body length of the resource fish. Equations (1.1–1.3) represent biomass dynamics and Eqs (1.4, 1.5) represent trait dynamics, where d⁢Sid⁢t=Vi⁢∂∂⁡Si⁢[d⁢Ci/d⁢tCi] and d⁢SRd⁢t=VR⁢∂∂⁡SR⁢[d⁢R/d⁢tR]. Here we assume that selection is frequency-dependent in both the resource (R) and consumers (C_i). τ_i is the maturation time for all individuals in a consumer population (McCann, 1998). Y_i is Beverton–Holt recruitment of age class i based on reproductive effort τ_i years ago and X_i is the growth of surviving recruits over τ_i, in which growth follows a Bertalanffy function (McCann, 1998). We assume the resource fish shows logistic growth: r is the intrinsic rate of increase of the resource, α is the intraspecific competition coefficient for the resource; Moreover, we use Holling type-II functional response to show predation terms: x_i and h_i are the attack rate and handling time of fish in age class i on the resource, respectively; d_i is the instantaneous rate of mortality of fish in age class i; δ_i is the instantaneous rate of reproductive energy invested into offspring for fish in age class i; a_i is the conversion costs of production of soma to gonadal tissue; g_i is the conversion efficiency of prey biomass into adult biomass; b_i is the egg supply rate of fish in age class i; d_Ji is density-independent egg mortality rate; d_JD is egg mortality related to egg density dependence; k_i is the rate of fish growth; W₁ is the asymptotic mass; and W₂ is the mass of an individual egg; both V_i and V₄ are the genetic variation of body length S_i and S_R, respectively; a higher value of the genetic variation means a faster speed of evolution. Finally, based on the relationship among resource and fish, similar to the early theoretical approach (Dercole et al., 2010), trait dependencies are modeled using the following functional forms:

Where w_i0 (i = 1, 2, 3) is the optimal trait values. r₀, r₁, x_i0, e_i, l_i, f_i, d_i0, and m_i are all positive. The growth rate r is maximum at S_R = w₄₀, where the growth of the resource is best adapted to its environment. The attack rate x_i is maximum at S_i = w_i0 and S_R = w₄₀ when the body length of cod matches with the prey fish. The mortality d_i of a fish in age class i is minimum at S_i = w_i0.

For eco-evolutionary model (Thompson, 2005), most of the parameter values were derived from the previous research (McCann, 1998) while the remaining set of parameters were obtained based on the reasonable guess values from parameter space presented in Supplementary Figure 3. All parameters can be found in Table 1. Similar to early theoretical work (McPeek, 2017; Cortez, 2018), we mainly explore how the speed of evolution in terms of the magnitude of genetic variation (V_i) affects the stability of the eco-evolutionary model (Thompson, 2005). We studied the dynamical properties of the eco-evolutionary model using numerical simulation method. Bifurcation diagrams and LE spectrum can be utilized to find the conditions under which the eco-evolutionary model produced chaotic dynamics (positive LEs). Moreover, we calculated both the ACF and LE of both empirical and theoretical time series to compare model predictions with field data. All simulations were carried out using MATLAB 7.0 (MathWorks, 2004).

The natural time series of body length suggests that chaotic dynamics may be presented in this fish species (Figure 1). First, we obtain the time delay τ = 2 (Figures 1D–F) and embedding dimension m = 11 (the rate of change of lnC(m, D) with lnD does not change with the increase of m; Figures 1G–I) for the age class I, II, and III, respectively. Then, by combining τ = 2 with m = 11, we obtain both of the LE values are greater than zero (Figures 1G–I). So, time series of the body length in the fish population exist in the chaotic region. Moreover, a similar result is obtained when the environmental noise of time series is filtered out (Supplementary Figure 2), all the LE values are larger than 0 (Supplementary Figures 2C,F,I), indicating dynamical chaos. So, the body length dynamics of the fish population present an intrinsic chaotic property.

As revealed by a theoretical simulation (Figure 2), the eco-evolutionary model can give rise to a rich set of potential dynamics. When the genetic variation of body length in the age class I, V₁, is small (V₁ = 0.03, LE = −0.0016; Figure 2A), the dynamics of body length S₁ approaches a simple, regular oscillation (period-2 oscillation). After an increase in the magnitude of genetic variation in body length S₁, V₁ (V₁ = 0.27, LE = −0.0012; Figure 2B), the dynamics of body length S₁ presents a doubling-periodic oscillation. Further, when V₁ is increased even more, the dynamics of body length S₁ appears to enter a chaotic regime (V₁ = 0.5, LE = 0.0643; Figure 2C) as suggested by the Lyapunov exponent (Figure 2F). Finally, periodic dynamics (regular oscillation) is observed when V₁ is increased even further (V₁ = 0.764, LE = −0.001; Figure 2D).

Since chaos may occur purely due to ecological dynamics, we did a bifurcation analysis assuming ecological dynamics only (i.e., assuming V₁ = V₂ = V₃ = V₄ = 0). Under these purely ecological conditions, the single-parameter bifurcation shows that the ecological process may not produce chaotic dynamics (Supplementary Figure 3), and the dynamics of population biomass C₁ is merely represented as a stable equilibrium.

Both theoretical predictions and empirical results show that the dynamics of a phenotypic trait, such as body length, can show chaotic properties. However, if the dynamics of our theoretical model is indeed a plausible representation of the empirical data, other dynamical properties of their respective time series should also be similar.

In this regard, the ACF analysis indeed indicates dynamical similarities among natural (Figure 3A) and theoretical (Figure 3B) time series. Both theoretical data and empirical data show two major peak lags: short (about 7 years) and long periods (about 24 years). However, due to the limited length of empirical data (in total 59 data points), the long period may be not significant (Figure 3A). Moreover, a predictable period (years) of the chaotic timeseries is 1/LE(S₃) (i.e., 1/0.1309 = 7.64 years in body length S₃; Figure 1I), which is close to 7 years (calculated by the ACF). Therefore, these results illustrate the eco-evolutionary model can depict the real dynamics of body length in a single fish system.

Early theoretical work found that regulating the magnitude of genetic variation in phenotypic trait can result in periodic oscillations (McPeek, 2017; Cortez, 2018) and chaotic dynamics (Dercole et al., 2010; Gilpin and Feldman, 2017). However, empirical evidence for eco-evolutionary chaos in natural ecosystems is lacking. In this study, we first use the phase space reconstruction method to present the evidence for potential chaos in fish body length. In our eco-evolutionary model, similar to recent theoretical approaches (Cortez, 2018; Yamamichi and Letten, 2021), we only analyze how system stability varies with the genetic variation (V₁). Our results show the increase of V₁ can generate chaotic dynamics in both phenotypic trait and biomass (LE > 0; Figure 2F). However, when we do not consider trait evolution (V_i = 0) in the model, population dynamics of the fish population may not present chaos. Then the ACF of both empirical data and simulation data of the eco-evolutionary model show that the predictable period (years) of body length in fish population is about 7 years. Therefore, the eco-evolutionary model can approximately depict and predict the intrinsic dynamic behavior of empirical data.

Moreover, recent theoretical research showed that ecological process (population density/biomass) will produce population oscillation if evolutionary process (phenotypic trait) has the oscillatory dynamics (McPeek, 2017; Cortez, 2018). In our study, we first present chaotic dynamics (irregular periodic oscillation) of the body length (attractor reconstruction), while comparative analyses of the ACF reveal broadly consistent results between experimental data and theoretical models. Therefore, we infer that the eco-evolutionary chaos may be driven by body length evolution (changes in genetic variation of body length) in the cod population.

In the eco-evolutionary model, the chaotic evolutionary trajectories (LE > 0; Figure 3F) are intrinsically unpredictable, and strongly dependent on the magnitude of genetic variation in phenotypic trait (V₁). The change switches between the non-chaotic (LE < 0) and chaotic dynamics (i.e., the “edge of chaos”) was discussed in earlier work (Ferriere and Fox, 1995; Turchin and Ellner, 2000; Benincà et al., 2015), which may balance the resource-fish system near the “evolutionary sliding” (Dercole et al., 2006). This maybe provide an evolutionary explanation for chaotic dynamics in the nature (Turchin, 2003). Our work identifies eco-evolutionary chaos in a natural ecological system. This empirical finding supports the theory that eco-evolutionary feedbacks can produce chaotic dynamics.

The rapid evolution of fish population has caused widespread discussion (Olsen et al., 2004; Jørgensen et al., 2007; Eikeset et al., 2016). In this study, we infer that the change of the body length (chaotic dynamics of phenotypic trait) is due to evolution. Moreover, the body length evolution has a limited predictability due to its chaotic property.

In the cod population, we checked for chaotic dynamics (positive LE) of body length based on both empirical data and eco-evolutionary model analyses. We found the LE of empirical data [LE(S₃) = 0.1309], which means predictable periods (1/LE) of empirical data (i.e., 7.64 years in body length S₃) is also close to 7 years (calculated by the ACF). The eco-evolutionary model can approximately depict and predict the intrinsic dynamic behavior of empirical data, that is, the limited predictable period of rapid evolution (V_i > 0) in body length of fish is approach to 7 years. Fish population growth is closely bound up with our daily life, considering the effectively predictable time of body length can help us make a suitable strategy for fishery conservation.