Limited effects of size-selective harvesting and harvesting-induced life-history changes on the temporal variability of biomass dynamics in complex food webs

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Limited effects of size-selective harvesting and harvesting-induced life-history changes on the temporal variability of biomass dynamics in complex food webs © 2022 The Authors. Published by Elsevier B.V. Published version Nonaka, Etsuko; Kuparinen, Anna Nonaka, E., & Kuparinen, A. (2023). Limited effects of size-selective harvesting and harvestinginduced life-history changes on the temporal variability of biomass dynamics in complex food webs. Ecological Modelling, 476, Article 110150. https://doi.org/10.1016/j.ecolmodel.2022.110150 2023 Ecological Modelling 476 (2023) 110150 Available online 8 December 2022 0304-3800/© 2022 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Limited effects of size-selective harvesting and harvesting-induced life-history changes on the temporal variability of biomass dynamics in complex food webs Etsuko Nonaka * , Anna Kuparinen Department of Biological and Environmental Sciences, University of Jyv¨ askyl¨ a, PO Box 35, 40014, Jyv¨ askyl¨ a, Finland ARTICLE INFO Keywords: Allometric trophic network (ATN) model Ecological networks Size-selective fishing Fisheries-induced evolution Niche model Temporal variation ABSTRACT Harvesting has been implicated in destabilizing the abundances of exploited populations. Because selective harvesting often targets large individuals, some studies have proposed that exploited populations often experience demographic shifts toward younger, smaller individuals and become more sensitive to environmental fluctuations. The theory of consumer–resource dynamics has been applied to address the impacts of harvesting in simple modular food webs, but harvested populations are embedded in a complex food web in nature. In addition, exploited populations have been shown to undergo trait evolution or phenotypic changes toward early maturation at smaller sizes. Using an empirically derived complex food web model, Kuparinen et al. (2016. Fishing-induced life-history changes degrade and destabilize harvested ecosystems. Scientific Reports: 6, 22,245) demonstrated that exploited fish populations and other ecosystem properties increased variability due to harvesting and harvesting-induced life-history evolution in the absence of environmental fluctuations. In the present study, we examined a large set of simulated complex food webs to attest the generality of the findings in Kuparinen et al. (2016). We found that harvesting both increased and decreased temporal variability of fish biomass dynamics, especially in food webs with intrinsically oscillating dynamics, while the vast majority of the food webs did not experience destabilization. We also elucidated the large impacts of the shape of functional responses on food web structure, energy flow, and changes in temporal variability caused by harvesting. Our results suggest that the destabilizing or stabilizing effects of harvesting and harvesting-induced evolution importantly depend on the shape of functional responses in dynamical models of complex food webs and that food webs with intrinsically oscillating dynamics are more prone to changes in temporal variability caused by harvesting and harvesting-induced evolution. 1. Introduction Size-selective harvesting, which preferentially removes individuals of certain sizes, has been implicated to result in phenotypic and potentially evolutionary changes in life-history traits related to survival and reproduction in exploited populations (Ernande et al., 2004; Heino et al., 2015; Kuparinen and Festa-Bianchet, 2016). It imposes high mortality and strong selection pressure typically on older and larger individuals and thus can truncate size structure (Barnett et al., 2017; Hsieh et al., 2006). The theory of life-history evolution predicts greater investment in reproduction and early maturation to evolve in populations undergoing high mortality in adults or equally across all life stages (resulting in smaller body size at maturity; Charlesworth, 1994; Law, 1979; Reznick and Ghalambor, 2005). As a consequence, populations experiencing high mortality in adults may consist more of younger and smaller individuals with higher reproductive capacities (Reznick et al., 1990; Stearns, 2000) and metabolic rate (Auer et al., 2018). The possibilities of fisheries-induced evolution (FIE) have been demonstrated by laboratory experiments of size-selective harvesting (Conover and Munch, 2002; Therkildsen et al., 2019; Uusi-Heikkil¨ a et al., 2017) and indicated by a diverse array of species in the wild (Allendorf and Hard, 2009; Darimont et al., 2009; Heino et al., 2015; Jorgensen et al., 2007; Kuparinen and Festa-Bianchet, 2016). Over the last several decades, reductions in the body sizes of fishes have been Abbreviations: ATN, Allometric trophic network; FIE, Fisheries-induced evolution. * Corresponding author. E-mail address: [email protected] (E. Nonaka). Contents lists available at ScienceDirect Ecological Modelling journal homepage: www.elsevier.com/locate/ecolmodel https://doi.org/10.1016/j.ecolmodel.2022.110150 Received 12 June 2022; Received in revised form 5 September 2022; Accepted 19 September 2022 Ecological Modelling 476 (2023) 110150 2 observed in exploited populations in response to size-selective fishing (Devine et al., 2012; McClenachan, 2009; Olsen et al., 2004). Evolutionary and plastic changes in life-history traits could modulate population growth rate, interaction strength, and energy flow up through food chains. Therefore, size-selective harvesting could alter the stability of populations (Bronikowski et al., 2002; Cameron et al., 2013) and of consumer–resource interactions (McCann, 2012; Rip and McCann, 2011). A modeling study by Kuparinen et al. (2016) showed that fishing-induced changes in life-history traits magnified temporal variation of intrinsic fish biomass dynamics and of other organisms over and above fishing even after fishing ceased, using an empirically parameterized bioenergetic food web model of a well-studied lake ecosystem. Although the notion that fishing often increases temporal variability of exploited fish populations appears to prevail (Anderson et al., 2008; Hsieh et al., 2006), the theory of consumer–resource interactions suggests that fishing could potentially either destabilize or stabilize fish populations through altered energy flow due to elevated metabolic rate and increased mortality of consumer populations (McCann, 2012; Rip and McCann, 2011). The generality of the effects of fishing and FIE on fish population dynamics suggested by existing studies remains unknown, albeit such information is highly relevant to sustainable management of fisheries ecosystems. Comprehensive empirical investigations are virtually impossible due to the lack of long-term fisheries and survey data that include not only measurements on relevant traits of focal species but also monitoring of the key prey species and the entire ecosystem dynamics. A complementary approach is to apply bioenergetic food web models, such as the one used by Kuparinen et al. (2016), to a large set of computer-generated (“synthetic”) realistic food webs that encompass a wide variety of network structures. This is the path we take in the present study. We investigated how typical it would be for fishing and fisheriesinduced life-history trait changes (hereafter FIE; these changes can be evolutionary or plastic responses) to increase temporal variabilities of food web dynamics and examined the associations between food web structure and responses to fishing and FIE in a large set of synthetic food webs. To this end, we used a recently developed composite model by Nonaka and Kuparinen (2021) which adds stage structure to fish populations in synthetic food web topologies stochastically generated by the niche model (Williams and Martinez, 2000). The dynamical model uses the allometric trophic network (ATN) approach (Bland et al., 2019; Boit et al., 2012; Brose et al., 2006; Glaum et al., 2020a; Hale et al., 2020; Martinez et al., 2006; Martinez, 2020; Williams et al., 2007; Williams and Martinez, 2004), which utilizes allometric scaling relationships for parameterizing vital rates, predator–prey interactions, and energy flow. We simulated deterministically without external forcing the scenarios for fishing with (F-FIE) and without FIE (fishing only; F-O) on an identical set of food webs along a gradient of the shape of the functional responses of consumers. Functional responses describe how the consumption rate of a predator varies across prey density and have critical effects in models of population dynamics. The parameter, the so-called Hill exponent (q in Eq. (2) below), shifts the shape of functional responses from type II (hyperbolic) to III (sigmoidal), along which the dynamical stability of consumer–resource interactions increase (Murdoch and Oaten, 1975; Uszko et al., 2015). Type III functional response promotes prey persistence at low prey densities through reduced predator consumption rates presumably owing to changes in prey and predator behaviors (e.g., prey hiding in refugia or predators forming search image on more abundant prey, adaptively foraging, or switching prey). It substantially influences the dynamical stability of ATN models (Martinez et al., 2006; Williams and Martinez, 2004). We classified food webs by types of responses to the F-FIE and F-O scenarios and measured food web structure and energy flow to compare food webs that showed different types of temporal dynamics and responses. 2. Methods 2.1. Life-history structured allometric trophic network model (ATN) The methods for generating food webs and the parameterization of the ATN are fully described in Nonaka and Kuparinen (2021; Fig. 1a). All the simulations were run in MATLAB with Statistical Toolbox and Parallel Computing Toolbox Release 2020b (The Mathworks, Inc., Natick, Massachusetts, United States). The code is available in the Dryad Digital Depository (https://datadryad.org/stash/dataset/doi:10.5061/dryad. z612jm6bk). Briefly, the model was built on the niche model (Williams and Martinez, 2000) with an extension to add life-history stages to a subset of “trophic species.” A trophic species consists of one or more taxa that share the same sets of predators and prey (Cohen et al., 1990; Havens, 1992). Some of the trophic species were grouped based on niche overlap to form a stage-structured population. The minimum overlap between feeding ranges to be qualified as consecutive stages of the same stage-structured fish taxon was set to 20% based on a literature survey (Nonaka and Kuparinen, 2021). Relative body masses are assigned to the terminal stages of stage-structured species and invertebrates based on the short-weighted trophic levels (Williams and Martinez, 2004) in accordance with the theory of allometric predator–prey body mass ratio (Brose et al., 2006). We used the von Bertalanffy isometric growth curve (Wootton, 1998) to define the masses of younger stage classes (Table 1). Because of this way of parameterization, the final predator–prey body mass ratios deviated from the theoretical mass ratio, but the shape of the distribution conformed to the empirical data (Brose et al., 2006, see the online appendix for Nonaka and Kuparinen 2021). It originally calculates body sizes (lengths) of fish as a function of age (Pauly 1980, Froese & Binohlan 2000), while we used ontogenetic stages instead because age and stages are practically related in many fishes (e.g., perch, white fish) and in our model (see below the use of the Leslie matrix). We applied the ATN to the food webs to simulate population biomass dynamics (Brose et al., 2006). The population dynamics of the species and life-history stages in the food web were described by a set of ordinary differential equations (ODE), which encoded the rates of biomass changes during a growing season. The population dynamics within the food webs were formulated as a multispecies consumer–resource model (Yodzis and Innes 1992; Brose et al., 2006; Boit et al., 2012; Bland et al., 2019), which was originally developed by Yodzis and Innes (1992) and has been further developed and analyzed by a number of studies (Bland et al., 2019; Boit et al., 2012; Brose et al., 2006; Glaum et al., 2020a; Hale et al., 2020; Martinez et al., 2006; Martinez, 2020; Williams et al., 2007; Williams and Martinez, 2004). Our parameterization followed the standard procedure for ATN models published elsewhere (Bland et al., 2019; Brose et al., 2006; Nonaka and Kuparinen, 2021). The population dynamics were encapsulated in the following equations. dBi dt =gi(1−∑ j∈autotrophs Bj K)Bi ⏞⏟⏟⏞ logistic growth of autotrophs −∑ j∈consumers xjyjiBj Fji eji ⏞⏟⏟⏞ loss to grazing (1) dBi dt = − fmxiBi ⏟⏞⏞⏟ metabolic loss +∑ j∈resources faxiyijBiFij ⏟⏞⏞⏟ dietary intake(′′gain′′) −∑ j∈consumers xjyjiBj Fji eji ⏟⏞⏞⏟ loss to predation −FmaxSstageBi ⏟⏞⏞⏟ loss to selective fishing . (2) Here, Bi was the biomass of species i, gi was the intrinsic growth rate of autotroph i, K was the carrying capacity, xi was the metabolic rate of consumer i, yij was the maximum consumption rate relative to metabolic rate, eij was the assimilation efficiency of predator i eating prey j, fm was the fraction of assimilated carbon lost for maintenance, and fa was the fraction of assimilated carbon that contributes to biomass growth (see Table 2a for parameter definitions and values). Fmax (maximum fishing intensity) and Sstage (stage-specific harvesting selectivity) together E. Nonaka and A. Kuparinen Ecological Modelling 476 (2023) 110150 3 determined the proportion of biomass removed by fishing (see2.2 Implementing fishing andfisheries-induced evolution). Fij was the functional response of consumer i when dealing with prey j (Bland et al., 2019; Boit et al., 2012): Fij = ω ij ∑l∈resources ω il Bq j Bq 0ij +∑k∈consumer(ckjpik BkBq 0kj )+∑l∈resources( ω il ∑l∈resources ω il Bq l)(3) ω ij was the diet preference parameter, and we assumed that fishes preferred fish to invertebrates and invertebrates to basal species as food items if they were included in their diets (Table 2a; Nonaka and Kuparinen 2021). ω ij was normalized by the sum of preferences toward all resources eaten by j to obtain relative preference. Likewise, we assumed that invertebrates preferred invertebrates over fish and basal species (Table 2a). Boij was the half saturation constant of consumer i eating prey j, which was the density of j at which iachieved half its maximum consumption rate (Table 2a; Williams, 2008). The Hill exponent q of the functional response was varied between 1.05 and 1.65 by an increment of 0.15. This range was chosen because q =1 and q =2 correspond to type II and III functional responses, respectively, and because past ATN studies used q between 1.2 and 1.5 (see 4.2. The effect of q (Hill exponent of functional responses) and the choices of q by previous studies; Williams and Martinez 2004b, Martinez et al., 2006). The second term in the denominator encoded interference competition between consumer iand kwhich shared the fraction pik of resource j at an intensity specified by the coefficient ckj (Table 2a). We assumed only intraspecific interference competition and no interspecific or between-stage interference competition for this study (i.e., ckj =0 for k∕= i; see 4.4. Limitation of the model in Discussion). Each year consisted of a 90-day growing season (90 time-steps in the ODE Solver) followed by a non-growing season where surplus energy (dietary intake minus all the losses in Eqn (2)) was allocated to growth and reproduction. The proportion allocated to reproduction increased linearly with stage and was none (0) for the first stage and 0.2 (Imax) for the terminal stage with an equal increment between stages (Table 2a; Kuparinen et al., 2016, Nonaka and Kuparinen 2021). Fish in the terminal stage additionally reproduced in exchange for 30% of their biomass (Roff, 1992; Wootton, 1998). The fraction of mature individuals of fish at each life stage was determined by a logistic function which specified that 50% of individuals were mature halfway to the terminal stage (Table 1). We used the Leslie matrix to shift the majority of somatic biomass to the next stage via growth and to convert it to new recruitment via reproduction at the terminal stage (Table 1). A fraction, ah, of a stage class remained in the same stage to reflect variation in growth (Table 2a). There was no gain or loss of biomass between consecutive growing seasons (i.e., no winter mortality). Fig. 1. a) A conceptual diagram of the model. This paper added selective fishing and fisheries-induced phenotypic changes (FIE) to the model developed by Nonaka and Kuparinen (2021) to analyze their impacts on fish food web dynamics. b) Simulation design. Persisting webs were the webs reached to the end of the pre-fishing period, all the nodes of persisting species/fish stages were connected in one network, no invertebrates were without predators, and at least one fish with 3 or more stages remained. F-FIE=fishing and fisheries-induced evolution. F-O =fishing only. E. Nonaka and A. Kuparinen Ecological Modelling 476 (2023) 110150 4 We generated 60,000 networks of 60 nodes (S0) and connectance (C) equal to 0.15, which is a commonly used value of connectance in ATN models (Berlow et al., 2009; Bland et al., 2019; Dunne et al., 2002; Martinez et al., 2012), with error tolerance at 0.025 (i.e., webs with C between 0.1475 and 0.1525 were accepted) containing between 2 and 6 fishes (Nfishes) with at least 3 and up to 5 stages (Nstagemin,Nstagemax; Table 2b). Taxa with trophic level (Thfish) >2 were eligible to become a fish node (i.e., need to include at least one non-basal species in diet). We then ran the ATN on each network from random initial biomasses uniformly distributed between 0.1 and 100 µgC/L. Taxa were considered extinct when the biomass was less than 10 −6 µgC/L or “invalid” if the biomass surpassed 10 12 µgC/L. This occurred because the terminal stage could accumulate large biomass over time (70% of biomass remained in the last stage) and some food webs allowed it to become very large. At the end of each growing season, fishes that retained only older stages, but not younger ones, for more than 10 years were removed as extinct. Food webs were regarded as persistent if simulations lasted for 1000 years, all the nodes of persisting species were connected in one network, no invertebrates were without predators, and at least one fish with 3 or more stages remained (“the qualifying criteria"). We refer to these selected webs as “persisting webs” below (Fig. 1b). 2.2. Implementing fishing and fisheries-induced evolution (F-FIE and F-O scenarios) On the persisting webs, we implemented two scenarios: (1) fishing and associated fisheries-induced life-history evolution (F-FIE) and (2) fishing only (F-O) (Fig. 1b). Following Kuparinen et al. (2016), we set stage-specific fishing selectivity (Sstage) to vary by stage, starting with zero in the first stage (Table 1). For example, for a fish species with 5 stage classes, Sstage is 0, 0.27, 0.73, 0.95, and 0.99 for stage-classes 1, 2, 3, 4, and 5, respectively. Non-fish trophic species were not subjected to fishing. Fishing effort (instantaneous fishing mortality rate, Fmax) was set at 0.5, as motivated by the empirical rate from Lake Constance used by Kuparinen et al. (2016). We imposed fisheries-induced changes in life histories, namely the onset of maturity and stage-specific metabolic rate, to investigate their effects on fish biomass dynamics. As in Kuparinen et al. (2016), during the fishing period, the age at 50% maturity in fishes was advanced by 0.5% per year (younger and smaller size at maturity) and fishes’ metabolic rate was increased by 10−4 per year (Andersen and Pedersen, 2009; Audzijonyte et al., 2013a; Darimont et al., 2009; Kuparinen et al., 2016). We assumed that evolutionary recovery from fishing was very slow (Conover et al., 2009; Enberg et al., 2009; Jorgensen et al., 2007; Kuparinen and Meril¨ a, 2007; Law, 2000; Law et al., 1989), so that fishes retained the life-history traits attained at the end of the fishing period unchanged during the post-fishing period. We first ensured that the biomasses of the species in the persisting webs settled at dynamic equilibria (stable or stable oscillations) and then subjected the fish populations to fishing and FIE for 50 years followed by 500 years of no fishing (F-FIE scenario). In the F-O scenario, no FIE was imposed during the fishing period (Fig. 1). The model was run for a long post-fishing period to observe the dynamics beyond the transients after fishing ceased. 2.3. Determining the temporal variability of fish biomass dynamics Coefficients of variation (CV) were calculated as the standard deviation of log10(biomass)divided by the absolute value of the mean across time for each fish stage. The maximum CV among stages was used to characterize the temporal variability of the fish biomass dynamics of a food web. We used the maximum, not the mean, to better capture the sensitivity of population fluctuations. CVs were calculated separately for the pre-fishing period from time step 901 to 1000 (CV pre ) and for the post-fishing period from time step 1451 to 1550 (CV post ) or from the last 50 time steps if the fish populations went extinct before the end of simulations. To address the main question regarding the additional effects of FIE more directly, we classified the simulated food webs by their responses to the F-FIE and F-O treatments. Because a set of identical food webs were subjected to the two scenarios, we considered their differential responses as isolated effects of FIE. We categorized the food webs into a 5 ×5 contingency matrix. (For the definitions of the response categories, see Table 3). The combinations of responses were further classified into “consistent” (same responses to F-FIE and F-O), “higher destabilization in F-FIE”, or “higher destabilization in F-O”. For the webs that were either “destabilized” or “stabilized” (hereafter refers to webs that increased and decreased temporal variability, respectively) consistently between the two scenarios, the levels of changes in temporal variability between F-FIE and F-O were quantitatively compared and assigned one of the three categories above (Table 3). Table 1 Model components. Model component Formulation Sources and notes Body mass at stage h (The von Bertalanffy isometric growth curve) Wh=W∞(1−e−κ(h−h0))3,where κ=3 v h∈ {1,…,v} v =terminal stage class of the fish (Wv W∞)=0.9 Pauly (1980), Froese & Binohlan (2000), Bland et al. (2019) The value of h0 is obtained by solving the equation for h0 with Wh= 0 and Wv from the predator-prey body mass ratio. The fraction of mature fish at each stage Pmature =1/(1+e−3(h−h50)) h50 =the stage at which 50% of individuals are mature Kuparinen et al. (2016) We assume h50 occurs halfway through to the terminal stage. Investment to reproduction I= (h− − 1)(Imax /(v−1)) Imax =maximum investment =0.2 Kuparinen et al. (2016) The Leslie matrix to model growth and reproduction by the terminal stage between growing seasons ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ Bi,1 Bi,2 Bi,3 ⋮ Bi,v ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠t+1 =⎛ ⎜ ⎜ ⎜ ⎜ ⎝ a10 0 0 bv b1a20 0 0 0b2a30 0 0 0 ⋱ ⋱ 0 0 0 0 bv−1av ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ Bi,0 Bi,1 Bi,2 ⋮ Bi,v ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠t bh =the proportion of biomass in stage h to be shifted to stage h+1 (or to stage 1 for h=v), ah =the proportion of biomass in stage hto remain in the same stage Modified from Bland et al. (2019) Younger reproducing stages reproduce through surplus energy, which is accounted for separately. Terminal stages reproduce through surplus energy and in exchange for biomass. Stage-specific harvesting selectivity Sstage =1/(1+e−2(h− ν /2)),for h>1 Sstage(h=1) = 0 Kuparinen et al. (2016) E. Nonaka and A. Kuparinen Ecological Modelling 476 (2023) 110150 5 Table 2 a. Parameter definitions and values used in this study for the allometric trophic network (ATN) model. Parameters Definition Value Sources or remarks Z Body mass ratio 102.6 for fish predators and prey 101.15 for invertebrate predators and prey Brose et al. (2006) gi Autotroph intrinsic growth rate for i Randomly drown from 0.8<N(0.9,0.5)<1 Bland et al. (2019) K Autotroph carrying capacity 540 μ gC/L Boit et al. (2012); Bland et al. (2019) xi Mass specific metabolic rate of i 0 for autotrophs 0.314M−0.15 i for invertebrates 0.88M−0.11 i for fish Bland et al. (2019); de Castro & Gaedke (2008); Killen et al. (2007, 2010) yij Maximum consumption rate of i eating j 4 for fish 8 for invertebrates Brose et al. (2006); Boit et al. (2012) eij Assimilation efficiency of i eating j 0.95 when jis a fish 0.75 when j is an invertebrate 0.45 when j is an autotroph Kelso (1972); Elliott (1976), Gavoni et al. (1986); Yodzis and Innes (1992); Brose et al. (2006) fm Fraction of assimilated carbon respired for maintenance of basic bodily functions 0.05 Modified from Bland et al. (2019); Boit et al. (2012) fa Fraction of assimilated carbon that contributes to growth 0.5 Modified from Bland et al. (2019); Boit et al. (2012) q Hill exponent of functional response (Eqn. (3)) Varied between 1.05 and 1.65 by an increment of 0.15 Williams & Martinez (2004); Martinez et al. (2006) ω ij Prey preference (relative to toward autotrophs) When iis a fish, 500 times toward fishes 200 times toward invertebrates When iis an invertebrate, 100 times toward invertebrates Fishes do not eat much autotrophs in temperate and northern regions (Gonzalez-Bergonzoni et al. 2012; Vejˇ ríkov´ a et al. 2016). Fish needs to eat fish to grow fast and large (Post 2003; Juanes et al. 2002). The values are calibrated to achieve little consumption of autotrophs by fishes with animal prey in the diets. Boij Half saturation density of ieating j 1.5 μ gC/L when an invertebrate i eats j 15 μ gC/L when fish i eats fish j 20 μ gC/L when fish i eats omnivore j 150 μ gC/L when fish i eats small herbivores j (more than 50 times smaller than the fish in body mass) 15 μ gC/L when fish i eats large herbivores j (not as much smaller than the fish in body mass) Tonin (2011); Martinez et al. (2012); Bland et al. (2019) See Fig. 1 in Bland et al. (2019) Herbivores are organisms whose diets consist of autotrophs for more than 70% (Bland et al., 2019). cij Consumer interference competition coefficient of i eating j Randomly drawn from 0 ≤Exp(λ=5) ≤ 0.5 when k is an invertebrate 3×10−4 when fish k eats fish j 10−4when fish k eats omnivore j 1 when fish i eats small herbivore j 10−4when fish k eats large herbivore j Tonin 2011; Martinez et al., 2012; Bland et al., 2019 See Fig. 1 in Bland et al. (2019) All interspecific interference competition for feeding on prey j, (ckj), was set to zero (i.e., ckj =0 for k∕= i; intraspecific (within-node) interference competition only). pij the fraction of resources of consumer i shared with consumer j (the number of prey i shares with j)/the number of i’s prey Bland et al. (2019) 1 −I Fraction of biomass invested to growth 1 for the first stage class 0.8 for the last stage class Kuparinen et al. (2016) I =investment to reproduction (Table 1) Pmature Probability of reaching maturation at stage h,h=1,2,…,n 0 for the first stage class (1+e−3(n−n50))−1 for higher stages, where n50 =n 2 is the stage at which 50% of the individuals become mature Kuparinen et al. (2016) ah Fraction of biomass staying in the same stage 0.05 for h=1,…,n−1 0.7 for h=v Fish in the terminal stage reproduces by converting 30% of total biomass to offspring regardless of surplus energy. bh Fraction of biomass moving to the next stage 1−ah i.e., b υ =0.3 L Number of days in a growing season 90 A typical length of a growing season in northern temperate regions (Kuparinen et al., 2016; Bland et al., 2019) Thextinct Extinction threshold 10−6 μ gC/L When biomass goes below this value, the population is considered extinct Thexplode Explosion threshold 1012 μ gC/L When biomass exceeds this value, the population is considered no longer sustainable. h50IncRate The rate of advances in stage at 50% maturity 5×10−3per year Kuparinen et al. (2016) and citations thereof. metIncRate The rate of increase in metabolism 10−4per year Kuparinen et al. (2016) and citations thereof. Fmax Maximum fishing intensity 0.5 (50% of fish biomass in target stages is harvested) Kuparinen et al. (2016) (continued on next page) E. Nonaka and A. Kuparinen Ecological Modelling 476 (2023) 110150 6 2.4. Quantifying food web structure and energy flow We used common network metrics in complex food web studies to quantify the basic structures of food webs (Table 4; Dunne, 2009; Gross et al., 2009; Williams and Martinez, 2008) and food web energetics by measuring the slopes of biomass spectra and the amounts and variation of energy influxes gained by fishes. Because we wanted to compare the food web structure before and after the fishing period (but not during the fishing period due to difficulties with removing transient dynamics), we measured these quantities with regard only to persisting species (i.e., “filtered” by the ATN dynamics) in the last 25 years of both the pre-fishing and the post-fishing periods (or the last 25 years before collapse). The slopes of biomass spectra were obtained by fitting a linear model with biomass in trophic levels as the response variable and body mass as the explanatory variable (Trebilco et al. 2013, Woodson et al. 2018). The model output the gain metrics containing the pairwise “gain” term (energy influxes) of the biomass rate equations Eqns. (1) and (2) integrated over a growing season. We took the mean across log10(gain)in the last 100 years of the pre-fishing and post-fishing periods (or minimum of 25 years before collapse) to calculate CVs, skewness (here skewness is invariant to scaling and shifting of distributions), and total sum across incoming links for all the fish stages. 3. Results 3.1. Temporal variability of fish biomass dynamics The proportion of persisting and qualified webs increased with q, while the proportion of oscillating qualified webs generally decreased (Fig. 2a). Food webs with stably oscillating biomass dynamics were rare. Almost the same number of food webs reached the end of simulations (i. e., 1550 years) in the F-FIE and F-O scenarios. The distributions of CVs of energy gain by fishes were highly skewed to the right, and the frequency rapidly and continuously declined toward larger CVs. Categorizing simulated time series with CV >0.05 as temporally variable (i.e., “unstable”) best matched visual inspection of subsets of the time series (Fig. S1). The distributions of absolute differences between CV pre and CV post were wide but mostly concentrated near 0 (Fig. S2). The absolute differences greater than 0.05 appeared the most Table 3 Classification of fish biomass dynamics in this study. Within each scenario (F-FIE or F-O) Stable (i.e., not temporally variable) CV <0.05 Unstable, oscillating CV >0.05 Destabilization by fishing (and FIE) (CV post – CV pre ) >0.05 Stabilization by fishing (and FIE) (CV pre – CV post ) >0.05 (De)stabilization after fishing or not No change in temporal variability abs(CV post – CV pre ) <0.05 Stabilization by fishing (and FIE) (CV pre – CV post ) >0.05 Destabilization by fishing (and FIE) (CV post – CV pre ) >0.05 Collapsed during or after fishing No fish with intact stage structure remained in these periods Contingency table (Fig. 3) Qualitative comparisons (No change, Collapse/during fishing, Collapse/after fishing) Consistent (blue) Same qualitative responses to F-FIE and F-O (in the diagonal) Higher destabilization in F- FIE (orange) F-FIE caused greater destabilization than F-O in the destabilization hierarchy* Higher destabilization in F-O (green) F-O caused greater destabilization than F-FIE in the destabilization hierarchy* Quantitative comparisons (Destabilized, Stabilized) Consistent (blue) abs(CV post , F-FIE – CV post,F-O ) <0.05 Higher destabilization in F- FIE (orange) (CV post , F-FIE – CV post,F-O ) >0.05 Higher destabilization in F-O (green) (CV post , F-O – CV post,F-FIE ) >0.05 * Destabilization hierarchy: Collapse/during >Collapse/after >Destabilized >No change >Stabilized. Table 4 Network metrics used in this study and their definitions. Name definition num_spp Number of trophic species (nodes) link_dens Link density (L/S =total number of links/total number of nodes) prop_fish_stages Proportion of fish stages prop_inverts Proportion of invertebrates prop_basal Proportion of basal species prop_intermed Proportion of intermediate species (species that have both predators and prey) prop_top Proportion of top species (species that do not have predators) prop_herbiv Proportion of herbivores (feed only on basal species) prop_cannib Proportion of cannibals (including inter-stage cannibals) prop_omniv Proportion of omnivores (species that feed on prey at different trophic levels) GenSD Standard deviation of generality (the number of prey) VulSD Standard deviation of vulnerability (the number of predators) max_TL Maximum trophic level ShPath_avg Average shortest path max_TL_fish Maximum trophic level of fish bm_spectra_slope Slope of biomass spectra fish_energy_influx Total amount of energy influxes into fishes gain_skew Skewness of the distribution of normalized energy influxes to fishes gain_CV Coefficient of variation of normalized energy influxes to fishes b. Parameter definitions and values used in the modified niche model to generate food webs (Nonaka and Kuparinen 2021). Parameter Definition Value used Sources or remarks S0 Number of species 60 C Connectance (proportion of realized links (L) out of all possible, L /S2 0) 0.15 Berlow et al. (2009); Martinez et al. (2012); Bland et al. (2019) Dunne et al. (2002) compiled empirically measured connectance and reported that it ranged between 0.03 and 0.3. Cerror Error tolerance on connectance 0.025 Bland et al. (2019) Nfishes The number of stage structured fish species Between 2 and 6 Reasonable numbers of naturally cooccurring fish species in a community Thfish A node at the trophic level >Thfish can become a fish stage 2 The diet of fishes should include non-autotrophs (trophic level of pure herbivores is 2) OLmin Minimum overlap of niche ranges between consecutive stages of a fish 0.2 In terms of the fraction of the union of the two feeding ranges Nstagemax The maximum number of stages a fish species can have 5 Reasonable maximum numbers of fish stages Nstagemin The minimum number of stages a fish species should have 3 Reasonable minimum numbers of fish stages Table 2 (continued) E. Nonaka and A. Kuparinen Ecological Modelling 476 (2023) 110150 7 adequate to indicate differences in temporal variability in comparisons with visual inspection (Fig. S3). Hence, we categorized food web dynamics as either “stable” (CV pre , CV post <0.05) or “unstable” (i.e., oscillating) (CV pre , CV post >0.05) and determined that fishing (and FIE for F-FIE) increased the variability of the dynamics if (CV post – CV pre ) > 0.05 and reduced if (CV pre – CV post ) >0.05 (Table 3). The majority of webs did not change temporal variability before and after the fishing period. The proportions of the webs that became more (destabilized) or less (stabilized) temporally variable were hence small but notably peaked at q=1.2 (Fig. 2b). The peak was dominated by food webs that were stable prior to fishing (Fig. 2b, c). Among the oscillating webs, the proportions of stabilized webs increased with the values of q up to 1.35, after which they stayed similar, while the proportions of destabilized webs remained roughly constant (Fig. 2d). The proportions of food webs that changed temporal variability, either increased or decreased, were greater under the F-FIE than the F-O scenario across all the values of q examined, indicating additional effects of FIE on the temporal variability of food web dynamics. Among the webs with oscillating dynamics, the proportions of the webs that were more negatively affected or more destabilized by F-FIE than by F-O were higher than those affected more by F-O than by F-FIE at lower values of q (Fig. 3). At higher values of q, the two became similar. For stable webs, most webs exhibited no change in both scenarios (Fig. S4). This result mirrored the above finding that FIE had additional effects on the temporal variability of food web dynamics beyond the effects of fishing alone at lower values of q. 3.2. Food web structure and energy flow The modified niche model generated food webs with a wide variety of network structural properties (Fig. S5). The values of the network structural metrics of the persisting food webs also showed high variation (Fig. 4). The mean trend almost monotonically increased or decreased along with q, and the patterns were similar among stabilized, destabilized, and unchanged webs in both the pre-fishing and post-fishing periods (Fig. 4a-o). As q increased, the persisting webs contained more species, more links, proportionally more intermediate and invertebrate species, and correspondingly proportionally fewer fishes and basal species. The patterns of the proportions of cannibals and omnivores and of maximum trophic levels reflected the compositional change (i.e., more fish stages). The destabilized and stabilized food webs on average had slightly lower proportions of basal species and slightly higher proportions of fishes that occupied marginally higher trophic levels than did the webs whose temporal variability not change before and after fishing (Fig. 4c,o). The mean trend for the slopes of biomass spectra increased with increasing q, indicating that the shapes of biomass pyramids became less bottom heavy as qincreased (Fig. 4p). The total energy gained by the fishes also increased with q (Fig. 4q). These increasing patterns together with increasing patterns of maximum trophic levels suggested that the food webs channeled more energy toward the top of the food chains (mostly fishes) with increasing q. The variation in the energy flow along each link into fishes (a proxy for interaction strength) increased with q (Fig. 4s), and an increasing number of these links became more concentrated at smaller values (greater positive skews; Fig. 4r) as q increased. These indicated that, when q was high, the food webs contained higher proportions of weaker feeding links on average. In the post-fishing period, food webs contained fewer species and links and experienced changes associated with these compositional changes (Fig. 4). The mean trends of the slopes of biomass spectra increased slightly (Fig. 4p) and the skewness and variation of energy fluxes into fishes declined modestly (Fig. 4r, s) at low to intermediate values of q for the stabilized and destabilized webs. Hence, fishing caused the food webs to have lower proportions of weak interactions and less bottom-heavy biomass pyramids at low to intermediate values of q. The patterns were very similar for the F-O scenario (Fig. S6). 4. Discussion 4.1. Limiting effects of fishing and FIE on the temporal variability of complex food webs Kuparinen et al. (2016) demonstrated that, using an empirically parameterized model without external forcing, increased temporal variability of fish abundances due to fisheries-induced phenotypic changes could remain even after fishing stopped. In the current study, Fig. 2. a) The proportions of qualified webs (i.e., reached 1000th year) containing at least one oscillating fish stage and of surviving webs (i.e., reached 1550th year), from the F-FIE and F-O scenarios. b) The proportions of destabilized (filled circle) and stabilized (open circle) webs among all the qualified webs. c) The proportions of destabilized webs that were stable in the pre-fishing period. d) The proportions of destabilized and stabilized webs that were oscillating in the pre-fishing period. E. Nonaka and A. Kuparinen Ecological Modelling 476 (2023) 110150 8 we examined their result in a large number of food webs with various topologies and observed that phenotypic changes caused by size-selective fishing infrequently increased or decreased but predominantly neither the temporal variability of fish biomass dynamics from the pre-fishing to the post-fishing period. The principle of energy flux (Gellner and McCann, 2016; McCann, 2012; Rip and McCann, 2011) can help explain the patterns of temporal variability in fish abundance. It predicts that an increase in energy flux to a consumer, relative to its loss term (relative energy flux), makes the population dynamics more temporally variable (or less stable) and the consumer–resource biomass ratio less bottom-heavy. Based on this principle, we can expect that, in the absence of environmental fluctuations, fishing and FIE could reduce the variability of fish populations if an increase in the loss term through fishing mortality and elevated metabolism reduce relative energy influx to the consumer. Likewise, we can expect that they could increase variability if energy influx to the consumer increases due to increased metabolic demands or higher resource availability mediated by trophic cascades or compensatory growth of prey. Finally, we can also expect that changes in gain and loss may be insufficient or that greater loss balances out increased gain, causing not much net changes in temporal variability. Our results appear to support the latter possibilities the most. In our study, the vast majority of food webs had stable biomass dynamics before and after fishing. Food web theory, which encompasses the principle of energy flux, that has been developed by McCann and others Gellner and McCann, 2016; McCann, 2012; McCauley et al., 2018; Nilsson et al., 2018; Rip and McCann, 2011) states that along the gradient of increasing energy flux (or decreasing energy loss) to a consumer a food web can transition from a non-excitable (never oscillates) to an excitable state (possible to oscillate), when the largest eigenvalue of the community matrix (the Jacobian matrix of the dynamical equations, Eqns. (1)-(3) turns from real to complex, before exhibiting sustained oscillations. In the stable food webs we observed, fishing and FIE presumably did not alter the food web structure or dynamics enough to bring them into the parameter region of sustained oscillations (Borrelli et al., 2015). Our result does not seem unrealistic as intrinsically oscillating populations, either exploited or not, appear uncommon in nature (Kendall et al., 1998; Shelton and Mangel, 2011). Instead, studies have pointed out the importance of environmental fluctuations interacting with nonlinear population dynamics perturbed by harvesting to destabilize the abundances of exploited populations (Anderson et al., 2008; Cameron et al., 2016). We observed that changes in temporal variabilities of fish biomasses occurred more often in food webs already oscillating prior to fishing and that those webs were more likely to reduce, rather than increase, temporal variability. The oscillating food webs were in excitable states (i.e., the largest eigenvalue of the community matrix was complex), where slight shifts along the energy flux gradient due to fishing and FIE could alter biomass variability qualitatively or quantitatively. The reason for reduced temporal variability could be that it is more likely for fishing and FIE to perturb the system to move toward a more stable state (by removing energy or pruning less stable modules), while it is less likely for perturbations to bring systems upward the instability gradient. Nonetheless, temporal variability increased in some oscillating webs because the systems collapsed or reorganized possibly to channel greater Fig. 3. a-e) The responses of identical food webs to the F-FIE and F-O scenarios for oscillating fish biomasses during the pre-fishing period were classified into the contingency matrices (for stable food webs see Fig. S4). The responses from F-FIE are listed along the vertical axis, and those from F-O are along the horizontal axis. The colored bars indicate the proportions of webs whose responses fell in the categories (the proportion goes up to 0.47 in all columns). The total number of food webs is shown at the bottom right corner. In the diagonal, F-FIE and F-O led to consistent responses in the first three categories (No change, Collapsed during fishing, Collapsed after fishing). In the Stabilized and Destabilized categories, the responses were further classified quantitatively with the same color scheme. f) The proportions of oscillating webs at five values of q on which FIE had an additional destabilizing effect (orange) and an additional stabilizing effect (green). E. Nonaka and A. Kuparinen