Limited predictability of body length in a fish population

Recent theoretical studies have identified chaotic dynamics in eco-evolutionary models. Yet, empirical evidence for eco-evolutionary chaos in natural ecosystems is lacking. In this study, we combine analyses of empirical data and an eco-evolutionary model to uncover chaotic dynamics of body length in a fish population (northeast Arctic cod: Gadus morhua). Consistent with chaotic attractors, the largest Lyapunov exponent (LE) of empirical data is positive, and approximately matches the LE of the model calculation, thus suggesting the potential for chaotic dynamics in this fish population. We also find that the autocorrelation function (ACF) of both empirical data and eco-evolutionary model shows a similar lag of approximately 7 years. Our combined analyses of natural time series and mathematical models suggest that chaotic dynamics of a phenotypic trait may be driven by trait evolution. This finding supports a growing theory that eco-evolutionary feedbacks can produce chaotic dynamics.


Introduction
One of the main questions in evolutionary biology is how ecological interactions affect the phenotypic trait evolution of species (Thompson, 2005). Variation in phenotypic traits (e.g., body size, behavior, morphology, and physiology) are a common feature in natural populations and phenotypic traits can evolve in natural environments (Coulson et al., 2011;Hanski, 2011;Agrawal et al., 2013). Mutual adaptation induced by ecological interactions; that is, coevolution, can shape the adaptive peaks of pairs of associated species (Guimarães et al., 2017). Moreover, selection caused by ecological interactions can even affect the dynamics of whole ecological systems (Koskella and Brockhurst, 2015) and accelerate the adaptation process of natural population (Galetti et al., 2013). These ecological and evolutionary interactions are critical to the structure and functioning of biodiversity (Ehrlich and Raven, 1964;Thompson, 2005). Hence, understanding how ecological interactions shape biodiversity will determine how coevolution acts on associated species (Iwao and Rausher, 1997;Parchman and Benkman, 2002;Ridenhour, 2005; Thompson et al., 2013).
The interplay of ecological and evolutionary dynamics has been an increasingly active area of research (Schoener, 2011;Koch et al., 2014). To this end, one of the most important areas of research is the long-term predictability of evolution (Green, 1991;Kauffman and Johnsen, 1991;Ferriere and Fox, 1995), where recent theoretical work has focused on the chaotic properties of ecological and evolutionary processes (Dercole et al., 2010;Schreiber et al., 2011;Gilpin and Feldman, 2017). For example, Dercole et al. (2010) coupled an evolutionary equation to a food chain model and found that coevolution can drive the population size to oscillate at the edge of chaos. Gilpin and Feldman (2017) showed that an even simple twospecies predator-prey model could give chaotic dynamics when evolution was included in the model. Schreiber et al. (2011) established an eco-evolutionary model showing how predator trait variation affects the trait value of prey, showing that rapid evolution of a polyphagous predator may generate chaotic dynamics due to eco-evolutionary feedbacks. Some researchers have also explored chaotic dynamics of multidimensional phenotypic traits and considered the effects of high-order interactions (e.g., the rate of phenotype x change with time (dx/dt) is affected by the x 3 term) on the rate of traits change with time, with the aim of analyzing the unpredictability of long-term evolution under high-order coupling (Doebeli and Ispolatov, 2014;Rego-Costa et al., 2018).
It is important to note that proving if eco-evolutionary feedbacks can lead to chaotic dynamics is a daunting task, and the empirical evidence for chaos in phenotypic traits are rare. This lack of evidence may be primarily due to the absence of long time series data on phenotypic traits. Interestingly, body size-a vital phenotypic trait-is also well known to be a major player in the dynamics and stability of interactions and here some longer time series exist (De Roos et al., 2003;Rooney et al., 2010;Heckmann et al., 2012;Delong et al., 2015). For example, Delong et al. (2015) analyzed the interaction between body size and trophic cascades, and found that the loss of larger predators have greater consequences on trophic control and biomass structure than smaller predators. Rooney et al. (2010) systematically analyzed the inner property of two food webs-the Cantabrian Sea Shelf marine and the Central Plains Experimental Range (CPER) shortgrass prairie soil. They found that biomass turnover rates (Production: Biomass ratio) decrease with increasing body size and larger organisms tend to have higher trophic positions. Heckmann et al. (2012) used a bioenergetics approach to analyze the interplay of body-size structure and adaptive foraging of consumers, and they found that stronger body-size structures (i.e., species on higher trophic levels have larger body masses than species on lower levels) and faster adaptation stabilize food webs. De Roos et al. (2003) synthesized research about population dynamics and body size dependence in individual life history, and they found that body size generally leads to population cycles driven by differences in competitiveness of differently sized individuals. Given that body size has been identified as important in mediating dynamical outcomes, it is surprising that no research has investigated ecoevolutionary dynamical properties (e.g., chaos) in the body size.
All of the previously mentioned research has shown that body size affects population dynamics and food web stability. Consistent with this, recent macroecological results found that body size from aquatic ecosystems (composed of smaller body-sized organisms) were less stable (a higher coefficient of variation) than wetland and terrestrial organisms (Rip and McCann, 2011). At the same time, for freshwater fish, recent work showed that the linear correlation between community biomass and mean body mass was not significant (Hatton et al., 2015). However, due to the influence of eco-evolutionary feedbacks (Ferriere and Legendre, 2013), we speculate that there is a non-linear correlation between community biomass and body mass. One well-known mechanism for eco-evolutionary feedbacks that may play a potent role in fish population dynamics is fisheries-induced evolution, which has been significant and ubiquitous in harvested ecosystems (Olsen et al., 2004;Jørgensen et al., 2007). Empirical data of a single fish population has found that body length of different ages can yield cyclic dynamics over time (Eikeset et al., 2016), it remains unknown whether these potentially and rapidly evolution induced evolutionary feedback can have chaotic properties or not.
In practice, the Lyapunov exponent (LE) is a classical method to characterize the chaotic nature of real and model ecosystems (Li et al., 2020;Rogers et al., 2022), and two approaches have been used to compute the LE: direct estimation and Jacobian/indirect estimation. On the one hand, the direct approach uses the definition to estimate LE from the data by measuring the divergence rate of the nearest neighbors over a finite time horizon (Rosenstein et al., 1993). In ecology, this approach is mainly used to describe the results of experimental systems (Becks et al., 2005;Graham et al., 2007;Benincà et al., 2008;Kosuta et al., 2008;Becks and Arndt, 2013;Wang et al., 2019). On the other hand, the Jacobian/indirect approach requires fitting a delayed embedding model (with embedding dimension and lag) to the available time series and calculating the LE from the Jacobian matrix of the model (Nychka et al., 1992). A variety of methods may be used to estimate unknown model frameworks, such as generalized additive models (GAMs) (Benincà et al., 2015), neural networks (Ellner and Turchin, 1995), local linear regression (Sugihara, 1994), and non-linear local LE (Li and Ding, 2022). In particular, to accurately predict the model framework of empirical data, the Jacobian/indirect method usually requires the incorporation of abiotic factors (e.g., temperature) into the system, which inevitably increases the data quality requirements.
In what follows, we use a biologically plausible ecoevolutionary model to study the relationship between a phenotypic trait (body length) and population biomass. We first examine dynamical properties of body length based on the direct LE estimation to reveal evidence of underlying oscillations [via autocorrelation function (ACF)] and chaotic dynamics (via attractor reconstruction) of the phenotypic trait in a natural time series. Second, we explore a relatively simple, yet plausible, theoretical eco-evolutionary model with body length as a trait, that shows similar dynamical signatures and chaos. Our results suggest that potential for eco-evolutionary dynamics. We answer two vital problems: how phenotypic trait affects population dynamics and system stability; and whether changes in the magnitude of genetic variation in body length can drive eco-evolutionary chaos (chaotic dynamics in ecological and evolutionary processes) in fish population dynamics. Our work supports a growing theory that eco-evolutionary feedbacks can produce chaotic dynamics.

Time series data
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.

Phase space reconstruction of body length
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.

Eco-evolutionary model
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(Cortez, , 2018McPeek, 2017;Yamamichi and Letten, 2021). The model can be written as: Frontiers in Ecology and Evolution 03 frontiersin.org Where i = 1, 2, 3, and (1.7) 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. 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 1 is the asymptotic mass; and W 2 is the mass of an individual egg; both V i and V 4 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 0 , r 1 , 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 40 , 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 40 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 .

Frontiers in Ecology and Evolution
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).

Empirical results
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.

Theoretical predictions of eco-evolutionary model
As revealed by a theoretical simulation (Figure 2), the ecoevolutionary model can give rise to a rich set of potential dynamics. When the genetic variation of body length in the age class I, V 1 , is small (V 1 = 0.03, LE = −0.0016; Figure 2A), the dynamics of body length S 1 approaches a simple, regular oscillation (period-2 oscillation). After an increase in the magnitude of genetic variation in body length S 1 , V 1 (V 1 = 0.27, LE = −0.0012; Figure 2B), the dynamics of body length S 1 presents a doubling-periodic oscillation. Further, when V 1 is increased even more, the dynamics of body length S 1 appears to enter a chaotic regime (V 1 = 0.5, LE = 0.0643; Figure 2C) as suggested by the Lyapunov exponent ( Figure 2F). Finally, periodic dynamics (regular oscillation) is observed when V 1 is increased even further (V 1 = 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 1 = V 2 = V 3 = V 4 = 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 1 is merely represented as a stable equilibrium.

Comparison between theoretical predictions and empirical data
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 3 ) (i.e., 1/0.1309 = 7.64 years in body length S 3 ; 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.

Discussion
The predictability of long-term evolution has puzzled ecologists for decades (Green, 1991;Kauffman and Johnsen, 1991;Ferriere and Fox, 1995). Although recent theoretical work has analyzed the chaotic properties of eco-evolutionary dynamics (Dercole et al., 2010;Schreiber et al., 2011;Doebeli and Ispolatov, 2014;Gilpin and Feldman, 2017;Rego-Costa et al., 2018), eco-evolutionary chaos has not been shown in empirical data. In this study, we combined an eco-evolutionary model with analyzes of empirical data to uncover chaotic dynamics of body length of a single fish population. Our work may prove the eco-evolutionary chaos in the natural ecological system.
Our work may be treated as a supplementary method to study fish species. In this study, we provided a new mechanism, that is, the genetic variation of phenotypic trait determines chaotic dynamics of body length. The new mechanism revealed Phase space reconstruction of single time series of body length. (A-C) Empirical data of body length; (D-F) Attractor reconstruction (τ = 2; based on the C-C method) of empirical data with noise; (G-I) Lyapunov exponent (LE) of body length: D is the neighborhood radius, m is the embedding dimension, and C is the correlation integral, which is influenced by both D and m. In each rectangle, m varies from 5 to 12 (from top line to bottom line, m = 11 can be obtained for each case when the rate of change of lnC(m, D) with lnD does not change with the increase of m). Body length S 1 (age class I: age 3-5), S 2 (age class II: age 6-9), and S 3 (age class III: age 10-12) are the mean of each age class. The original data is shown in Supplementary Figure 1. in our work is different from recent studies (Eikeset et al., 2016;Andersen, 2019), which showed many mechanisms could influence population dynamics (stable or oscillatory state) of fish species, such as the change of age structure, lifehistory parameters, and fishing. Actually, our work involved with above factors. First of all, for different study methods (attractor reconstruction and eco-evolutionary model), different age classes (age 3-12) were divided into three broader age classes, and each class has the different maximal attack rate and reproductive rate (Table 1). Moreover, life-history parameters (i.e., individual growth rate k i and reproduction rate b i in Table 1) indeed influence population dynamics (non-chaos in Supplementary Figure 3) of fish species when we neglected the evolutionary dynamics (V 1 = V 2 = V 3 = V 4 = 0). Finally, although we do not introduce human harvest (fishing) into the eco-evolutionary model, we infer that fishing will impact the demography and recruitment of a fish stock. Therefore, it is worth trying to analyze how each ecological parameter influence population dynamics of fish. The simulation result showed that non-chaos emerges when we do not consider phenotypic trait evolution in the mechanistic model (Supplementary Figure 3).
In short, in this work, the free-equation approach (attractor reconstruction) proved chaos in trait dynamics, and then a plausible interpretation mechanism (regulating the genetic variation of phenotypic trait may result in chaotic dynamics of body length in the eco-evolutionary model) for chaos in body length of fish was suggested, that is, rapid evolution in phenotypic trait may result in chaos in both ecological (population density/biomass) and evolutionary (phenotypic trait) processes.

The road to eco-evolutionary chaos
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 1 ). Our results show the increase of V 1 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 1 ). 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.

Rapid evolution in body length of fish is limited predictability
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 Chaotic dynamics is observed for the northeast Arctic cod (Gadus morhua) based on empirical data and simulation data analyzes. (A) Autocorrelation function (ACF) of empirical data; (B) ACF of simulation data; body length class I (S 1 ), class II (S 2 ), and class III (S 3 ). Blue line denotes significance threshold and the value in the middle of two blue lines can be considered not significantly different from 0. 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 3 ) = 0.1309], which means predictable periods (1/LE) of empirical data (i.e., 7.64 years in body length S 3 ) 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.

Data availability statement
The original contributions presented in this study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Author contributions
LW designed the study, finished model simulations, carried out data analyses, and wrote the first draft of the manuscript. Both authors contributed substantially to revisions and approved the submitted version.

Funding
This research was supported by the National Natural Science Foundation of China (32201259) and Doctor Start-up Fund Project of Jiangxi Normal University (12021105).