Unusually Paced Life History Strategies of Marine Megafauna Drive Atypical Sensitivities to Environmental Variability

Understanding why different life history strategies respond differently to changes in environmental variability is necessary to be able to predict eco-evolutionary population responses to change. Marine megafauna display unusual combinations of life history traits. For example, rays, sharks and turtles are all long-lived, characteristic of slow life histories. However, turtles also have very high reproduction rates and juvenile mortality, characteristic of fast life histories. Sharks and rays, in contrast, produce a few live-born young, which have low mortality rates, characteristic of slow life histories. This raises the question if marine megafaunal responses to environmental variability follow conventional life history patterns, including the pattern that fast life histories are more sensitive to environmental autocorrelation than slow life histories. To answer this question, we used a functional trait approach to quantify for different species of mobulid rays, cheloniid sea turtles and carcharhinid sharks – all inhabitants or visitors of (human-dominated) coastalscapes – how their life history, average size and log stochastic population growth rate, log(λs), respond to changes in environmental autocorrelation and in the frequency of favorable environmental conditions. The faster life histories were more sensitive to temporal frequency of favourable environmental conditions, but both faster and slower life histories were equally sensitive, although of opposite sign, to environmental autocorrelation. These patterns are atypical, likely following from the unusual life history traits that the megafauna display, as responses were linked to variation in mortality, growth and reproduction rates. Our findings signify the importance of understanding how life history traits and population responses to environmental change are linked. Such understanding is a basis for accurate predictions of marine megafauna population responses to environmental perturbations like (over)fishing, and to shifts in the autocorrelation of environmental variables, ultimately contributing toward bending the curve on marine biodiversity loss.


INTRODUCTION
Being able to accurately predict how populations of organisms respond to environmental change is one of the key challenges for biologists today (Clements and Ozgul, 2016;Salguero-Gómez et al., 2016). How populations respond to changes in their environment is mediated by the survival, growth and reproduction rates of individuals Silvertown, 1996, 2004;Tuljapurkar and Haridas, 2006). Different combinations of these demographic rates comprise different life history strategies (Gaillard et al., 2016;Salguero-Gómez et al., 2016;Paniw et al., 2018). Different life histories, in turn, are linked to different types of population responses to environmental change. For example, slow life histories are characterized by low mortality, low fecundity and low development rates (Gaillard et al., 1989;Salguero-Gómez et al., 2016), and populations of slow life history are generally buffered against increased environmental variability (Morris et al., 2008(Morris et al., , 2011Dalgleish et al., 2010;Saether et al., 2013; but see Jongejans et al., 2010;McDonald et al., 2017). Fast life histories are characterized by high mortality, high fecundity and high development rates (Gaillard et al., 1989;Salguero-Gómez et al., 2016). In contrast to slow life history populations, fast life history populations are very sensitive to environmental variability, and their population sizes can fluctuate greatly if environmental variability is high (Morris et al., 2008;McDonald et al., 2017). Unraveling how life history strategies and demographic rates are linked can therefore provide in depth understanding of how populations respond to environmental change (Salguero-Gómez et al., 2016;Paniw et al., 2018;Smallegange and Berg, 2019).
Responses to environmental change of populations characterized by different life history strategies are extensively studied in stochastic population analyses, in which demographic rates and population growth vary over time (Getz and Haight, 1989;Lande et al., 1997). In these analyses, environmental autocorrelation (see Table 1, which has definitions of terms used in this study) is assumed to be absent; that is, the future environmental state is unrelated to the current one. In nature, however, environmental autocorrelation is usually positive; that is, the current and future environmental state are likely to be similar (Ariño and Pimm, 1995;Halley, 1996;Inchausti and Halley, 2002). Yet, we are only beginning to unravel whether links between individual survival, growth and reproduction rates and life history strategies identified in the absence of environmental autocorrelation hold when environmental autocorrelation is positive (or negative; Paniw et al., 2018). Environmental autocorrelation can leave a signature in the autocorrelation of demographic rates (Ariño and Pimm, 1995;Paniw et al., 2018). Theoretical studies have shown that a change in the autocorrelation of demographic rates may increase or decrease population growth rate depending on the life history strategy (Tuljapurkar et al., 2009). Using a cross-taxonomical approach, Paniw et al. (2018) illustrated that fast life histories were more sensitive to simulated autocorrelation in demographic rates than slow life histories. Smallegange and Berg (2019), on the other hand, found, in a species-specific analysis, that

Term Definition
Demographic rates A demographic rate expresses the rate at which an individual moves through its life cycle and includes an individual's rate of survival, growth and reproduction.

Environmental autocorrelation
The temporal correlation in environmental conditions.
Environmental noise High, positive values of environmental autocorrelation are referred to as red noise; high and negative values as blue noise; and when environmental correlation equals zero, this denotes white noise where the probability of switching states is independent of the current state.
Environmental stochasticity (Un)predictable, temporal fluctuation in environmental conditions. (Un)predictability is defined as an (in)ability to predict the future environmental state from the current state.

Environmental variability
The variability of environmental conditions over time.
Greater environmental variability is often associated with increased variability of resource supply.

Good environment frequency
The temporal frequency with which good environmental conditions occur.
Life history strategy The relative effects of demographic rates on how population growth responds to environmental change are captured in combinations of life history traits and comprise a population's life history strategy.

Stochastic environment
A stochastic environment is random in nature and not deterministic.
the slow life history reef manta ray Mobula alfredi was more sensitive to simulated environmental autocorrelation than the fast life history beach hopper Orchestia gammarellus. Not only environmental autocorrelation determines how demographic rates and population growth vary over time, also good environment frequency (Table 1): the temporal frequency with which favorable environmental conditions occur (Caswell, 2001). For example, Smallegange and Berg (2019) found that the fast life history species O. gammarellus was very sensitive to good environment frequency, in contrast to the slow life history species M. alfredi. With predicted global changes in environmental patterning (García-Carreras and Reuman, 2011) due to, e.g., shifts in environmental autocorrelation or good environment frequency, it is urgent to gain in-depth understanding of the life history processes that result in distinct demographic responses between different life history strategies to such shifts. Marine megafauna are globally recognized as providing significant economic, cultural, and ecological values (Bakker et al., 2016;Teh et al., 2018). Despite this, many marine megafauna species have experienced population declines, putting some species at risk of extinction (Paleczny et al., 2010;Davidson et al., 2012;Randhawa et al., 2015). Charismatic marine megafauna that inhabit coastalscapes are mobulid rays, carcharhinid sharks and cheloniid sea turtles. They fulfill key roles as predators or grazers in coastal ecosystems in the (sub)tropics (Ferretti et al., 2010;Roff et al., 2016), yet their persistence is threatened by (over)fishing, loss of breeding habitat, pollution, habitat disturbance and pathogens (Clarke et al., 2006;Estes et al., 2011;Brooks et al., 2013;Vianna et al., 2013;Ward-Paige et al., 2013;Dwyer et al., 2020;Jatmiko and Catur Nugroho, 2020). As a result, most of these marine megafauna are now listed as vulnerable, near-threatened or endangered on the international union for conservation of nature red list of threatened species (Marshall et al., 2011;Brooks et al., 2013).
Here we use a life history trait approach to explore if the subtly different life history strategies of mobulid rays, carcharhinid sharks and cheloniid sea turtles are differentially sensitive to environmental autocorrelation and good environment frequency. We describe individual life histories from a dynamic energy budget (DEB) perspective, and examine how the environment affects the change in individual life histories, thereby generating the dynamics of population structure (Webb et al., 2010). Such an approach can offer mechanistic insights into how individual level processes affect population responses to environmental change. The traditional approach to do so is to use physiologically structured population models (PSPMs). However, the analysis of PSPMs requires rather complicated methodology and PSPMs typically are more mathematical representations of biological systems (deRoos and Persson, 2013). We therefore resort to using the recently developed DEB integral projection model (DEB-IPM; Smallegange et al., 2017). In contrast to PSPMs, IPMs are data-driven, provide a way of synthesizing complex life history information and can be analyzed using more straightforward mathematical techniques (Smallegange and Coulson, 2013). The DEB-IPM takes data on individual life histories as input parameters to describe growth and reproduction following a simple version of the standard model of Kooijman's DEB theory (Kooijman, 2010), also known as the Kooijman-Metz model (Kooijman and Metz, 1984). Mortality and the association between parent and offspring characteristics do not follow DEB theory, and are estimated from individual-level observations (Smallegange et al., 2017).
Using the DEB-IPM, we examine population-level responses in fitness and in mean size (as a proxy for population structure) to a wide range of environmental autocorrelation and good environment frequency values, and ask (i) if the overall slower mobulid and carcharhinid life histories show higher sensitivities of log stochastic population growth rates, log(λ s ), to environmental autocorrelation than the turtles that also display some fast life history characteristics, (ii) if the different life histories show different sensitivities of log(λ s ) to good environment frequency, and (iii) if the log(λ s ) responses can be captured by any of the life history traits within the DEB-IPM (mortality, growth, reproduction) along an axis of life history speed. We parameterized DEB-IPMs for four species of mobulid rays, five species of carcharhinid sharks, and four species of cheloniid sea turtle. For each species, we carried out stochastic simulations to assess if they show high or low sensitivity of log(λ s ) and of mean body size to shifts in environmental autocorrelation and good environment frequency. Environments with varying environmental autocorrelation and good environment frequency were simulated using a stochastic demographic model in which the temporal sequence of good and bad food environments is driven by a Markovian process that governs the serial correlation of environment states. In order to cover a wide range of such stochastic environments ( Table 1), we varied this serial correlation from blue to white to red environmental noise color (corresponding to a negative first-order autocorrelation, no autocorrelation, and a positive first-order autocorrelation of the temporal sequence, respectively; Table 1) for two values of good environment frequency. Next, we varied good environment frequency from almost zero to almost one, in the absence of environmental autocorrelation (white noise: Table 1). Finally, across stochastic environments, we conducted perturbation analyses to assess which functional trait had the strongest effect on log(λ s ), and we assessed how functional traits relate to variation in log(λ s ). Our approach allows us to assess whether the effects of environmental autocorrelation and good environment frequency on a select group of marine megafauna demography can be predicted from the life history traits that comprise each species' life history strategy.

Brief Description of the DEB-IPM
The demographic functions that describe growth and reproduction in the DEB-IPM are derived from the Kooijman-Metz model (Kooijman and Metz, 1984), which is a simple version of the standard model of Kooijman's DEB theory, but still fulfils the criteria for general explanatory models for the energetics of individuals (Sousa et al., 2010). The Kooijman-Metz model assumes that individual organisms are isomorphic (body surface area and volume are proportional to squared and cubed length, respectively). The rate at which an individual ingests food, I, is assumed to be proportional to the maximum ingestion rate I max , the current feeding level Y and body surface area, and hence to the squared length of an organism: I = I max YL 2 . Ingested food is assimilated with a constant efficiency ε. A constant fraction κ of assimilated energy is allocated to soma (metabolic maintenance and growth); this energy equals κεI max YL 2 and is used to first cover maintenance costs, which are proportional to body volume following ξ L 3 (ξ is the proportionality constant relating maintenance energy requirements to cubed length), while the remainder is allocated to somatic growth. The remaining fraction 1-κ of assimilated energy, the reproduction energy, is allocated to reproduction in case of adults and to the development of reproductive organs in case of juveniles, and equals (1 − κ)εI max YL 2 . This means that, if an individual survives from year t to year t + 1, it grows from length L to length L following a von Bertalanffy growth curve, dL dt =ṙ B (L m · Y-L), where r B is the von Bertalanffy growth rate and L m = κεI max /ξ is the maximum length under conditions of unlimited resource. Both κ and I max are assumed to be constant across experienced feeding levels, and therefore L m is also assumed constant. If a surviving female is an adult, she also produces offspring. According to the Kooijman-Metz model, reproduction, i.e., the number of offspring produced by an individual of length L between time t and t + 1, equals Y · R m · L 2 /L 2 m . The parameter R m is the maximum reproduction rate of an individual of maximum length L m . Note that R m is proportional to (1-κ) (Kooijman and Metz, 1984), whereas L m is proportional to κ, which controls energy conservation. However, the role of κ in the DEB-IPM is mostly implicit, as κ is used as input parameter only in the starvation condition (see below), whereas R m and L m are measured directly from data. Like L m , R m is also proportional to I max ; since both κ and I max are assumed to be constant across experienced feeding levels, R m is also assumed constant.
The above individual life history events are captured in the DEB-IPM by four fundamental functions to describe the dynamics of a population comprising cohorts of females of different sizes (Smallegange et al., 2017): (1) the survival function, S (L(t)) (unit: y −1 ), describing the probability of surviving from year t to year t + 1; (2) the growth function, G(L , L(t)) (unit: y −1 ), describing the probability that an individual of body length L at year t grows to length L at t + 1, conditional on survival; (3) the reproduction function, R (L (t)) (unit: # y −1 ), giving the number of offspring produced between year t and t + 1 by an individual of length L at year t; and (4) the parent-offspring function, D(L , L(t)) (unit: y −1 ), the latter which describes the association between the body length of the parent L and offspring length L (i.e., to what extent does offspring size depend on parental size). The DEB-IPM assumes no effect of temperature on fundamental functions. Denoting the number of females at year t by N (L, t) means that the dynamics of the body length number distribution from year t to t + 1 can be written as: where the closed interval denotes the length domain. Implicitly underlying the population-level model of eqn 1, like in any IPM, is a stochastic, individual-based model, in which individuals follow Markovian growth trajectories that depend on an individual's current state (Easterling et al., 2000). This individual variability is in standard IPMs modeled in the functions describing growth, G(L , L(t)), and the parent-offspring association, D(L , L(t)), using a probability density distribution, typically Gaussian (Easterling et al., 2000). In the DEB-IPM, this individual variability arises from how individuals experience the environment; specifically, the experienced feeding level Y follows a Gaussian distribution with mean E(Y) and standard deviation σ(Y). It means that individuals within a cohort of length L do not necessarily experience the same feeding level due to demographic stochasticity (e.g., individuals, independently of each other, have good or bad luck in their feeding experience).
The survival function S (L(t)) in Eq. 1 is the probability that an individual of length L survives from year t to t + 1: where E(Y) can range from zero (empty gut) to one (full gut).
Here, individuals that experience E(Y) = 1 can be assumed to always have a full gut. Individuals die from starvation at a body length at which maintenance requirements exceed the total amount of assimilated energy, which occurs when L > L m · E(Y)/κ and hence, then, S (L (t)) = 0 (e.g., an individual of size L m will die of starvation if E(Y) < κ). Juveniles and adults often have different mortality rates, and, thus, juveniles (L b ≤ L < L p ) that do not die of starvation (i.e., L ≤ L m · E(Y)/κ) have a mortality rate of µ j , and adults (L p ≤ L ≤ L m ) that do not die of starvation (i.e., L ≤ L m · E(Y)/κ) have a mortality rate of µ a .
The function G(L , L(t)) is the probability that an individual of body length L at year t grows to length L at t + 1, conditional on survival, and, following common practice (Easterling et al., 2000;Coulson, 2012;Merow et al., 2014), follows a Gaussian distribution: with the growth realized by a cohort of individuals with length L(t) equaling (4) assuming individuals do not shrink under starvation conditions (Smallegange et al., 2017), and the variance in length at time t + 1 for a cohort of individuals of length L as where σ(Y) is the standard deviation of the expected feeding level.
The reproduction function R (L(t)) gives the number of offspring produced between year t and t + 1 by an individual of length L at year t: Individuals are mature when they reach puberty at body length L p and only surviving adults reproduce (Eq. 1); thus, only individuals within a cohort of length L p ≤ L ≤ L m Y/κ reproduce. When maintenance costs cannot be covered (starvation), ingested energy is rechannelled from reproduction to maintenance, which occurs for L > L m E(Y).
The probability density function D L , L(t) gives the probability that the offspring of an individual of body length L are of length L at year t + 1, and hence describes the association between parent and offspring character values: is the expected size of offspring produced by a cohort of individuals with length L(t), and σ 2 L b (L (t)) the associated variance. For simplicity, we set E L b (L (t)) constant and assumed its associated variance, σ 2 L b (L (t)), to be very small.

Parameterization of the DEB-IPM
We took values for the fraction of assimilated energy allocated to maintenance and growth, κ, from the Add My Pet database (Add-my-pet, 2020), if available, and otherwise assumed κ = 0.8, as in the generalized animal (Kooijman, 2010; Table 2); mortality rates are all estimates of natural mortality rates. Most parameter estimates were taken directly from the literature (references in Table 2), but for some species they had to be calculated from other life history parameters, which is explained below.

Mobulid Rays
In case of M. alfredi, we calculated juvenile mortality rate µ i from its yearling survival rate (P y = 0.63 yr −1 ), its juvenile survival rate (P j = 0.95 yr −1 ) and its average age at maturity (α = 10; Kashiwagi, 2014;Smallegange et al., 2016): Table 2), assuming P = e −µ (Caswell, 2001). In case of M. birostris, we were unable to find estimates for µ i . We instead took M. alfredi yearling survival rate (P y = 0.63 yr −1 ) and juvenile survival rate (P j = 0.95 yr −1 ; Kashiwagi, 2014;Smallegange et al., 2016), and, given the fact that M. birostris matures on average at 9 years of age (Dulvy et al., 2014;Rambahiniarison et al., 2018), µ i = −log 9 P y × P 8 j ( Table 2). We also could not find an estimate for M. birostris R m and therefore assumed it equaled that of M. alfredi (Table 2). Finally, in case of M. thurstoni, we were unable to find estimates for µ j and µ a , and assumed these to equal those of M. japanica ( Table 2).

Carcharhinids Sharks
In case of C. limbatus, N. brevirostris, and G. cuvier we calculated juvenile mortality rate, µ j , as µ j = −log α √ P m with P m as the survival probability from newborn to maturity, and α as average age at maturity (yr). In case of C. limbatus, P m = 0.26 and α = 7 yr; in case of N. brevirostris, P m = 0.12 and α = 12.7 yr, and in case of G. cuvier, P m = 0.02 and α = 9 yr (Ward-Paige et al., 2010; Table 2). For the latter three species, adult mortality rate, µ a , was calculated as µ a = −log l √ P a , where P a is the natural mortality rate over the adult life span, l. For C. limbatus, P a = 0.25 and l = 18 yr; for N. brevirostris, P a = 0.18 and l = 25 yr, and for G. cuvier, P a = 0.68 and l = 28 yr ( Table 2). In case of T. obesus and C. amblyrhynchos, µ j and µ a were calculated from life tables of age-specific, yearly survival rates (Robbins, 2006). Table 2).

Cheloniid Turtles
In case of C. mydas, C. caretta and E. imbricata we calculated µ i as µ j = −log α √ P m . In case of C. mydas, P m = 0.01 [as 10 out of 1,000 eggs need to survive for a population to persist (Hirth and Schaffer, 1974)] and α = 35 yr (Chaloupka and Limpus, 2005); in case of C. caretta, P m = 0.025 [as 2.5 out of 1,000 eggs need to survive for a population to persist (Frazer, 1986)] and α = 26 yr (Casale et al., 2011); and in case of E. imbricata, P m = 0.01 (we assumed, like C. mydas, that 10 out of 1,000 eggs need to survive for a population to persist) and α = 30 yr (Limpus, 1992;Chaloupka and Limpus, 1997; Table 2). For all cheloniid species, adult mortality rate, µ a , was calculated as µ a = −log (P s ), where P s is the annual adult survival rate. For D. coriacea, P s = 0.89 (Dutton et al., 2005); for C. mydas, P s = 0.95

Stochastic Demographic Model
We used the stochastic demographic model where p(t) is the population vector at time t and A (t) is a DEB-IPM at time t defined by a two-state Markov chain that gives the probability distribution of environment states at time t. In this chain, state 1 is the good environment and state 2 is the bad environment. This results in the following Markov chain habitat transition matrix H (Caswell, 2001, p. 379): where p is the probability of switching from the good to the bad environment, and q is the probability of switching from the bad to the good environment. The serial or autocorrelation of the Markov chain equals ρ = 1-p-q (Caswell, 2001, p. 379). High, positive values of ρ are referred to as red noise; high and negative values of ρ as blue noise; and ρ = 0 denotes white noise where the probability of switching states is independent of the current state ( Table 1). The good environment frequency is given by f = q/(p + q) (Caswell, 2001, p. 379). We used a high feeding level E(Y) high and a low feeding level E(Y) low to define the good and bad environmental states, respectively. Specifically, E(Y) low is the expected feeding level associated with population decline (λ < 1) and E(Y) high as the expected feeding level associated with population increase (λ > 1). For the rays and sharks, we set E(Y) low = 0.7 and E(Y) high = 1.0 and for the turtles that have a slightly faster life history speed, we set E(Y) low = 0.6 and E(Y) high = 1.0 (Smallegange and Berg, 2019). We set σ(Y) = 0.3. We aimed to independently study the effects of environmental autocorrelation ρ and good environment frequency f ( Table 1). We thus varied ρ across the full noise gradient, while keeping f fixed at f = 0.5; and over a gradient of white and red noise, while keeping f fixed at f = 0.75. We varied f across the full gradient ranging from almost zero to almost unity, while keeping ρ constant at ρ = 0.5. Each stochastic simulation was generated in MatLab (MATLAB. version 9.2.0.556344 [R2017a], The MathWorks, Inc., Natick, MA, United States) by iterating H over a time series of length 50,000 (with an initial transient length of 500 discarded, a starting population of one individual in each size bin, and with the initial environment state chosen randomly [see also Tuljapurkar et al., 2003]; e.g., Smallegange et al., 2014). This sequence determines the environment state, and hence the feeding level E(Y), that a population experiences at each time step, from which the individual-level functions were calculated to construct A (t), i.e., the DEB-IPM at time t defined by the feeding level E(Y) at time t. At each time step, A (t) was stored with associated vectors of population structure to calculate the log of the stochastic population growth rate λ s as log(λ s ) = 1 τ τ−1 τ=0 r t with r t = log p(t + 1)/p(t) , where τ = 50,000 -500 = 49,500. At each time step, also the mean of the body size distribution was calculated, after which we used the pooled mean body size, calculated as the grand mean of all mean body sizes over time period τ, for our analysis. For a given species, a decrease in pooled mean body size across a stochastic gradient would indicate an increase in the proportion of juveniles in the population, where an increase in pooled mean body size would indicate an increase in the proportion of adults in the population.

Perturbation Analysis
We conducted a perturbation analysis in MatLab (MATLAB. version 9.2.0.556344 [R2017a]) to examine the elasticity of the log stochastic population growth rate, log(λ s ), to perturbation of each of the life history parameters: length at birth L b , length at puberty L p , maximum length L m , juvenile mortality rate µ j , adult mortality rate µ a , von Bertalanffy growth rateṙ B , and maximum reproduction rate R m ( Table 2). This analysis allows us to identify which life history parameter is most influential to log(λ s ), and if this depends on the type of stochastic environment. In order to do this, we perturbed each parameter by 1% and calculated the elasticity of log(λ s ) to each model parameter. We excluded the parameter κ(fraction of assimilated energy allocated to maintenance and growth), because it cannot be perturbed directly as, apart from occurring in the starvation condition (Eq. 2), it is implicitly included in the model within L m (which is mathematically proportional to κ) and R m (which is mathematically proportional to [1-κ]; Kooijman and Metz, 1984).

Linking Life History Parameters to log(λ s ) Across Stochastic Environments
Four of the model parameters are life history characteristics that place species on the fast-slow life history continuum. Specifically, fast (or slow) life histories are characterized by high (or low) mortality rate µ j and µ a , high (or low) individual growth ratė r B , high (or low) reproduction rate R m (Gaillard et al., 1989: Stearns, 1992Salguero-Gómez et al., 2016). It would thus be interesting to assess if the magnitude change in log(λ s ) over the environmental autocorrelation, ρ, or good environment frequency, f, gradient is linked to species-specific parameter values. For example, the reef manta ray M. alfredi, a slow life history species, characterized by low values of µ j , µ a ,ṙ B , and R m , is more sensitive to a change in environmental autocorrelation than the fast life history species O. gammarellus, characterized by high values of µ j , µ a ,ṙ B , and R m (Smallegange and Berg, 2019). In a cross-species comparison, you would then expect a correlation between the values of µ j , µ a ,ṙ B , and R m and the magnitude with which the log stochastic population growth rate, log(λ s ), responds to changes in ρ. To this end, we first normalized (scaled)ṙ B and R m to compare across species (µ j , µ a already only take values between zero and unity), and then used a linear model with each (normalized) model parameter as a continuous explanatory variable to test its relationship to the difference between log(λ s ) at the end and at the start of each stochastic gradient. Analyses were carried out in MatLab (MATLAB. version 9.2.0.556344 [R2017a]).

Shifts in log(λ s ) and Pooled Mean Body Size Across Stochastic Environments
Across the gradient of environmental autocorrelation ρ, the log stochastic population growth rate, log(λ s ), of mobulids increased as ρ increased from zero (white noise) to high positive (red noise), both when good environment frequency f is fixed at f = 0.5 and f = 0.75 ( Figure 1A). For all mobulid species, log(λ s ) was higher when f is fixed at f = 0.75 (gray lines in Figure 1A) then when f = 0.5 (black lines in Figure 1A). Carcharhinids log(λ s ) values were hardly affected by environmental autocorrelation, except for the tiger shark (G. cuvier), for which log(λ s ) showed a humpshaped response over the environmental autocorrelation gradient ( Figure 1B). Like the mobulids, log(λ s ) of the carcharhinids was higher when f is fixed at f = 0.75 (gray lines in Figure 1B) than when f = 0.5 (black lines in Figure 1B), particularly for the tiger shark G. cuvier. Cheloniid log(λ s ) was the least sensitive to temporal autocorrelation, because log(λ s ) for most species only decreased very slightly over the environmental autocorrelation ρ gradient before showing a small but sudden drop in value as ρ approached unity ( Figure 1C). Like the mobulids and most carcharhinids, cheloniid log(λ s ) was higher when good environment frequency f is fixed at f = 0.75 (gray lines in Figure 1C) than when f = 0.5 (black lines in Figure 1C). Across the gradient of good environment frequency f (while keeping ρ = 0), log(λ s ) increased for all species, but the increase was greater for carcharhinids and cheloniids, than for mobulids (Figures 1D,E).
Pooled mean body size of mobulids and carcharhinids (Figures 2A,B) was largely insensitive to shifts in environmental autocorrelation ρ, both when good environment frequency f is fixed at f = 0.5 and f = 0.75 (Figures 2A,B), and was also largely insensitive to shifts in good environment frequency f (Figures 2D,E). The exception to the latter patterns was the tiger shark (G. cuvier), for which pooled mean body size increased with increasing ρ, both when f is fixed at f = 0.5 and f = 0.75 (Figure 2B), and decreased with increasing f (Figure 2E). In contrast, pooled mean body size of cheloniids increased with increasing values of environmental autocorrelation ρ, with a sudden, sharp increase as ρ approached unity, both when f is fixed at f = 0.5 and f = 0.75 (Figure 2C), whereas it decreased with increasing values of f ( Figure 2F). Finally, overall, pooled mean body size across the environmental autocorrelation ρ gradient was for some species, like the tiger shark (G. cuvier) and leatherback turtle (D. coriacea) higher when f = 0.75, than when f = 0.5, although the extent of this difference varied between species (Figures 2A-C).

Elasticity of log(λ s ) to Perturbation of Life History Parameters Across Stochastic Environments
The perturbation analyses revealed which life history parameter elicited relatively the highest change in the log stochastic population growth rate, log(λ s ), within a single stochastic environment. Four parameters influenced log(λ s ) most: an increase (or decrease) in length at puberty, L p , decreased (or increased) log(λ s ); an increase (or decrease) in maximum length, L m , increased (or decreased) log(λ s ); an increase (or decrease) in juvenile mortality rate, µ j , decreased (or increased) log(λ s ); and an increase (or decrease) in the von Bertalanffy growth rate,ṙ B , increased (or decreased) log(λ s ).
Along the environmental autocorrelation gradient of ρ, while keeping good environment frequency, f, fixed at f = 0.5, log(λ s ) was predominantly most sensitive to perturbation of maximum length, L m , across the mobulid species (Table 3). Exceptions were log(λ s ) of the reef manta ray M. alfredi, which was most sensitive to perturbation of length at puberty L p , and log(λ s ) of the bentfin devil ray M. thurstoni, which was most sensitive to perturbation of L p at the most negative values of ρ (Table 3). Carcharhinid log(λ s ) was, across species, most sensitive to perturbation of either L p or L m across the ρ gradient (Table 3). Cheloniid log(λ s ) was across most of the ρ gradient most sensitive to perturbation of L p and L m , but also to perturbation of juvenile mortality rate µ i and von Bertalanffy growth rateṙ B , depending on species and ρ value ( Table 3).
Running the same perturbation analyses across the environmental autocorrelation gradient ρ, but with the good environment frequency f fixed at f = 0.75, showed that the sensitivity response of mobulid log(λ s ) was very similar to when f = 0.5: mobulids were still mostly sensitive to perturbation of maximum length L m ( Table 3). Exceptions were again M. alfredi and M. thurstoni, for which, at a larger range of negative values of ρ, log(λ s ) is most sensitive to perturbation of length at puberty, L p . Carcharhinid log(λ s ) was equally sensitive to perturbation of L p or L m when f = 0.75 compared to when f = 0.5 (Table 3). Finally, cheloniid log(λ s ) showed the largest difference in sensitivity response across the environmental autocorrelation ρ gradient, depending on whether f = 0.5 or f = 0.75. When f = 0.75, none of the cheloniids log(λ s ) were sensitive to perturbation of juvenile mortality rate µ i , or the von Bertalanffy growth rateṙ B ; instead, log(λ s ) was always sensitive to perturbation of length at puberty L p , except for the loggerhead turtle, C. caretta, which log(λ s ) was always most sensitive to perturbation of maximum length L m ( Table 3).   Shown are the life history parameter to which log(λ s ) was most sensitive when perturbed by 0.01% (see main text). Parameters are: length at puberty L p (cm); maximum length L m (cm); juvenile mortality rateµ i (yr −1 ); and von Bertalanffy growth rateṙ B (yr −1 ).
Running the perturbation analyses across the gradient of good environment frequency, f (while keeping environmental autocorrelation ρ = 0) revealed that log(λ s ) showed strikingly similar patterns in elasticity to perturbation of life history parameters as observed across the ρ gradient, with f kept constant at f = 0.5 (Table 3). Specifically, each mobulid log(λ s ) was most sensitive to perturbation of the same life history parameter as observed across the environmental autocorrelation ρ gradient (with f = 0.5; Table 3). The higher values of M. thurstoni log(λ s ) were also mostly sensitive to perturbation of maximum length, L m , at higher values of f (log(λ s ) was also most sensitive to perturbation of maximum length, L m , at higher ρ values, with f fixed at f = 0.5; Table 3). For each carcharhinid species, log(λ s ) was mostly sensitive to perturbation of the same life history parameter as observed across the environmental autocorrelation ρ gradient (with f = 0.5; Table 3), with some small differences. At very low good environment frequency f values, log(λ s ) of the blacktip shark C. limbatus and whitetip reef shark T. obesus were most sensitive to perturbation of juvenile mortality rate, µ i , instead of length at puberty, L p , at the most negative ρ values (Table 3). Finally, cheloniid log(λ s ) was again most sensitive to perturbation of L p , L m or µ i , but notṙ B (as across the ρ gradient with f = 0.5: Table 3), as loggerhead turtle C. caretta and hawksbill turtle E. imbricata log(λ s ) were most sensitive to perturbation of µ i at low values of f (as opposed toṙ B along the lowest range of the ρ gradient with f = 0.5; Table 3).

DISCUSSION
Our goal was to assess whether the effects of environmental autocorrelation and good environment frequency on a select group of marine megafauna demography can be predicted from the life history traits that comprise each species' life history strategy. We found that responses in how the log stochastic population growth rate, log(λ s ), varied across the environmental autocorrelation and good environment frequency gradients were grossly captured by variation in three life history traits: juvenile mortality rate, (scaled) maximum reproduction rate and von Bertalanffy growth rate. Higher values of these traits are characteristic of faster life histories, whereas lower values of these traits are characteristic of slower life histories (Gaillard et al., 2016;Salguero-Gómez et al., 2016). Along this mortalityreproduction axis of life history speed we found that the turtles, the blacktip shark and the tiger shark displayed a higher increase in log(λ s ) with increasing good environment frequency than the other sharks and the rays. Across the environmental autocorrelation gradient from blue to red noise, we found that log(λ s ) of the turtles and the tiger shark decreased, whereas log(λ s ) of the other sharks and the rays increased in similar magnitude across the same gradient.

Are There General Life History Patterns Across Stochastic Environments?
Changes in the patterning of environmental stochasticity ( Table 1) may have significant implications for population viability worldwide (Heino and Sabadell, 2003;Ruokolainen et al., 2009;Fey and Wieczynski, 2016). It is thus urgent to gain an indepth understanding of the life history processes that result in distinct population responses between slow and fast life histories to shifts in environmental stochasticity. The conventional approach to do so partitions variation in life history statistics used to describe both ecological and evolutionary dynamics -like generation time, age at maturity and mean lifetime reproductive success -across many species along an axis of life history speed and an axis of reproductive strategy (Stearns, 1989, Gaillard et al., 2016Salguero-Gómez et al., 2016;Paniw et al., 2018;Capdevila et al., 2020). A population's position in this 2D life history space is then taken to inform on its response to environmental change. This conventional approach is phenomenological and lacks a mechanistic representation of the biological processes that give rise to the observed variation in life history statistics (Musick, 1999;Salguero-Gómez, 2017;Salguero-Gómez et al., 2018). Instead, we used a recently developed functional trait approach where the demographic processes of growth and reproduction are described mechanistically by individual energy budgets (Smallegange et al., 2017). Because the model is constructed from (functional) life history traits, we were able to directly relate variation in life history traits to the demographic consequences of shifts in environmental stochasticity.
Unsurprisingly, we found that, as the good environment frequency increased, the log stochastic population growth rate, log(λ s ), of all species increased, although species at the slow end of our mortality-reproduction axis of life history speed (predominantly the rays) showed a lower increase than the species at the fast end of the life history speed axis (particularly the turtles, blacktip shark, and tiger shark). We surmise that this difference in response is predominantly due to differences in reproduction rates between the slower and faster life history strategies. The faster life histories have relatively high maximum reproduction rates (R m ). Because for non-starving adults, reproduction is proportional to the product of expected feeding level E(Y) and R m (Eq. 6), an increase in the frequency of high E(Y) has a proportionally larger effect on population growth of fast life histories than that of slow ones. This effect is reflected in the response of pooled, mean body size across the gradient of good environment frequency. Whereas pooled, mean body size of the slower life histories showed at most a slight, steady decline across this gradient, pooled, mean body size of the faster life histories decreased greatly with increasing good environment frequency. The latter shift in population sizestructure reflects a substantial, relative increase in the proportion of small (juvenile) individuals compared to large (adult) ones, and thus an increase in reproduction, in line with the relatively strong increase in log(λ s ) over this gradient. Despite the fact that pooled, mean body size decreased with increasing good environment frequency, it should be noted that the largest size that an individual can attain, L ∞ (Kooijman, 2010), will increase because it is positively related with feeding level following the relationship L ∞ = L m .E(Y) (Smallegange et al., 2017).
A more surprising finding was that, over the environmental autocorrelation gradient from blue (high, negative values of ρ) to red noise (high, positive values of ρ), log(λ s ) of the faster life histories decreased in value, whereas log(λ s ) of the slower life histories increased in value of similar magnitude. This result is surprising because it contrasts with the result that Paniw et al. (2018) found, in which log(λ s ) of faster life histories changed in absolute value more over the environmental autocorrelation gradient from white to red noise than log(λ s ) of slower species. The latter pattern is in line with empirical and theoretical findings (Franco and Silvertown, 2004;Morris et al., 2008Morris et al., , 2011Salguero-Gómez et al., 2016), whereas the pattern that we found is not. One reason could be that the taxa that we studied differ in their life history traits in ways that do not fit the typical pace-oflife axis of variation in life history strategy that has emerged from large(r), cross-taxonomical studies on animals and plants (Gaillard et al., 1989;Salguero-Gómez et al., 2016;Paniw et al., 2018), including aquatic species (Cortés, 2002;Frisk et al., 2005;Quetglas et al., 2016;Capdevila et al., 2020). Specifically, the fast-slow life history strategy continuum is bounded by the fastliving end where species develop quickly, have high reproduction rates, but also high mortality rates, and the slow-living end where species develop slowly, have low reproduction rate, but also low mortality rates (Stearns, 1992). Following this characterization, empirical and theoretical findings have shown that slow life history species are typically less sensitive to environmental autocorrelation than fast life history species, precisely because they are long-lived (Franco and Silvertown, 2004;Salguero-Gómez et al., 2016), have long generation times (Tuljapurkar et al., 2009), and low juvenile and adult mortality rates (Morris et al., 2008(Morris et al., , 2011. The slower life history species in our study indeed have lower mortality rates than the faster species, but all the species that we studied are long-lived with long generation times. It is perhaps this atypical combination of life history traits that marine megafauna display that underlies the contrasting demographic responses to environmental autocorrelation of the slower and faster life histories, which, in our study, are characterized by low, respectively, high juvenile mortality and reproduction. The next step is thus to mechanistically understand how temporal stochasticity in environmental conditions drives these contrasting demographic responses between the turtles and tiger shark on the one hand, and the rays and other sharks the other hand. In response to a shift from blue to red environmental autocorrelation, log(λ s ) of all rays showed the highest increase in value, whereas log(λ s ) of the turtles and the tiger shark showed the highest decrease in value. The turtles and the tiger shark have in common that they have very high juvenile mortality rates, higher than those of the rays. Additionally, the turtles show by far the greatest leap in growth from hatchling to length at puberty. All of this means that any prolonged period of unfavorable environmental conditions, which occurred either half the time (f = 0.5) or a quarter of the time (f = 0.75) under the red noise environments we investigated, have high negative impact on juvenile persistence in these species. Long periods of unfavorable environmental conditions slow individual growth (Eqn 6), prolonging the period over which juveniles are exposed to high mortality until they reach their size at puberty. As a result, fewer will survive as adults to reproduce, lowering log(λ s ). This process is reflected in the population size-structure: pooled, mean body size of the turtles and tiger shark increased over the environmental autocorrelation gradient as log(λ s ) decreased, reflecting reduced reproduction and proportional decrease of smaller sized individuals. The rays, in contrast, appeared to profit from the prolonged periods of favorable conditions under red noise, judging from the increase in log(λ s ) and corresponding slight decline in pooled, mean body size (due to increased reproduction) across the environmental autocorrelation gradient from blue to red. These results provide a fresh perspective onto the theoretical prediction that populations that recover slowly from past perturbations [like long-lived [iteroparous] species in constant environments (Salguero-Gómez et al., 2016)] should be more sensitive to environmental autocorrelation than those that are more resilient (Tuljapurkar and Haridas, 2006). Our finding that long-lived marine megafauna are sensitive to environmental autocorrelation, coupled with the fact that a similar functional trait analysis revealed that the manta ray M. alfredi was more sensitive to environmental autocorrelation than the fast life history amphipod O. gammarellus (Smallegange and Berg, 2019), are in line with the latter theoretical prediction. Conventional approaches, however, have found that population persistence of long-lived animal and plant species in autocorrelated environments is buffered from environmental variation, and, crucially, this patterns was independent of reproductive strategy (e.g., Metcalf and Koons, 2007;Morris et al., 2008). Yet, we find that it is precisely the reproductive strategy of marine megafauna that plays such an important role in their demographic responses, because it determines whether populations grow or decline over the environmental autocorrelation gradient in environmental conditions. The finding by Paniw et al. (2018) that fast life histories were on average more sensitive to environmental autocorrelation than slow species explained 50% of the variance in sensitivity to environmental autocorrelation. Perhaps the remaining variation is unexplained because not all species (groups) strictly map onto the fast-slow life history continuum (this study; Jongejans et al., 2010;McDonald et al., 2017). Finding general life history patterns in how species respond to environmental change is important as it aids our prediction of population persistence, extinction, and diversification (Salguero-Gómez et al., 2016), particularly in the absence of species-specific information in case of many marine organisms (Heppell et al., 1999). However, we should not lose sight of the fact that general life history patterns are not always expected (Galipaud and Kokko, in press), in which case it can be crucial to mechanistically understand how specific interactions between environmental conditions and life history traits drive a population's response to environmental autocorrelation.

Some Implications for Assessing the Impact of Environmental Change on Marine Megafauna Persistence
The global demand for marine animal products such as shark fins (Clarke et al., 2006), swim bladders (Sadovy and Cheung, 2003;Clarke, 2004), and ray gill plates (White et al., 2006;Ward-Paige et al., 2013) is unsustainable (Berkes et al., 2006;Lenzen et al., 2012). Particularly for the slower life history species has the intense fishing exploitation that targets these demands resulted in population declines and increased risks of extinction, sometimes with synergistic effects of environmental conditions (Jennings et al., 1999;Schindler et al., 2002). Assessing how responses to environmental change of species with contrasting life history strategies differ within an ecological community is thus essential to manage mixed fisheries (Musick, 1999;Jennings and Rice, 2011;Link, 2013). Scholars of chondrichthyanscartilaginous fishes including rays and sharks -often use the maximum intrinsic rate of population increase, r max , to assess a population's status. When population trajectories are lacking, which is often the case for oceanic species (Bradshaw et al., 2007), r max can be a useful statistic to evaluate a species' relative risk of overexploitation by fishing (Dulvy et al., 2014) as it can be taken as the equivalent of the fishing mortality that drives a species to extinction (Myers and Mertz, 1998). The statistic r max is estimated by solving the Euler-Lotka equation (Myers and Mertz, 1998), the most simple versions of which take survival and reproduction schedules as input parameters (Stearns, 1992). Like the Euler-Lotka equation (Stearns, 1992), a DEB-IPM is parameterized by survival and reproduction rates, but also takes the von Bertalanffy growth rate and three length parameters as input (Smallegange et al., 2017). However, whereas the Euler-Lotka equation assumes density-independence (Stearns, 1992), a DEB-IPM can account for density-dependence, either simulated via the expected feeding level E(Y), or incorporated mechanistically by including resource dynamics (Smallegange et al., 2017). Many applications of the Euler-Lotka equation to chondrichthyan demographic responses to environmental perturbations set survival to maturity close to unity (García et al., 2008;Hutchings et al., 2012;Dulvy et al., 2014, but see Pardo et al., 2016a, assuming very high juvenile survival rates because chondrichthyans invest highly into offspring. Yet even within this group of species, reproductive strategies can differ greatly (Branstetter, 1990). Our analysis using the DEB-IPM warns against assuming high juvenile survival across chondrichthyans, because the different reproductive strategies of at least some carcharhinids sharks and mobulid rays can result in markedly different demographic responses to environmental change.
Our perturbation analyses revealed that the log stochastic population growth rate, log(λ s ), of all marine megafauna that we studied was, across most, if not all of the length of each stochastic gradient, most sensitive to perturbation of length at puberty, L p , or maximum length, L m . Because we took a life history approach, we can examine how a change in either parameter affects population performance. A decrease in L p , all else being equal, means that individuals start reproducing at a smaller size (Kooijman and Metz, 1984), but they may also produce fewer or smaller offspring (Hume, 2019). In rays and sharks, a decrease in size at maturity can reduce fecundity, because smaller mothers produce smaller offspring (Sibly et al., 2018;Hume, 2019), which have been postulated to have lower fitness (Motta et al., 2007;Hume, 2019). However, in, e.g., the lemon shark, smaller juveniles have higher survival rates than larger individuals of the same age, which could favor maturation at a smaller size (Dibattista et al., 2007). Additionally, several conspecific sea turtle populations differ in length at puberty (Goshe et al., 2010;Bell and Pike, 2012;Snover et al., 2013;Avens et al., 2015Avens et al., , 2017Avens et al., , 2020Marn et al., 2019). To investigate how the interaction between the benefits and costs of maturing at a smaller size affects population growth rates, we re-ran all our analyses for a scenario in which length at puberty of all cheloniid species is increased by 10% or decreased by 10% (Supplementary Appendix). We found that our results are robust against perturbation of length at puberty, because increasing or decreasing length at puberty of sea turtles by 10% did not qualitatively affect how different life history traits are linked to population responses to shifts in environmental stochasticity (Supplementary Appendix). The other parameter that log(λ s ) was very sensitive to was maximum length, L m . A decrease in L m will result in smaller individuals within a population (Eqn 4). Smaller individuals lay fewer eggs (because R m is related to L m (Smallegange et al., 2017: Supplementary Appendix) so that a decrease in L m can reduce population size. In the marine environment, a decrease in L m can be caused by the selective (over)fishing of large individuals, because prolonged (over)fishing of the largest individuals can impose selection on the developmental processes underlying growth and development and drive contemporary evolutionary responses in L p and L m toward earlier maturation at smaller sizes (Waples and Audzijonyte, 2016). In many fishes, selective (over)fishing has reduced mean body size (Frisk et al., 2005;Fenberg and Roy, 2008). For example, mean body size of whale sharks Rhincodon typus has declined in response to fishing, and, at the same time, population abundance has reduced (Bradshaw et al., 2009). These observations, in concert with our findings, signify the importance of understanding the eco-evolutionary interaction between (evolutionary) shifts in life history traits and population growth to be able to better predict population responses of marine megafauna to fishing and environmental change.

CONCLUSION AND OUTLOOK
We found that, across a range of mobulid, carcharhinid and cheloniid species, faster life histories were more sensitive to temporal frequency of good environment conditions, but both faster and slower life histories were equally sensitive, although of opposite sign, to environmental autocorrelation. These patterns are atypical, likely following from the unusual life history traits that these megafauna display. Our analysis is a first exploration of marine megafaunal life history strategies across stochastic environments. As such, we did not take into account specifics of the different species life histories, although our post hoc perturbation of length at puberty of sea turtles did not qualitatively affect how marine megafaunal life history traits are linked to shifts in environmental stochasticity (Supplementary Appendix). Most sea turtles undergo a major ontogenetic habitat shift between oceanic and neritic foraging areas (Ramirez et al., 2015(Ramirez et al., , 2017Tomaszewicz et al., 2017Tomaszewicz et al., , 2018. The disconnect and timing of occupation of different feeding habitats could have an impact on sea turtle population responses to shifts in environmental stochasticity. Observed changes in marine systems due to climate change include shifts in range and changes in algal, plankton and fish abundances (Intergovernmental Panel on Climate Change [IPCC], 2007), thus affecting food availability (Hamann et al., 2013) and potentially the environmental autocorrelation of food availability over time, or the frequency with which food availability is favorable (good environment frequency). Because such climate-induced changes can impact population responses to harvesting (Isomaa et al., 2014;Smallegange and Ens, 2018), future work should focus on understanding the impact of shifting environmental stochasticity on (harvested) populations (Huntingford et al., 2013;Boulton and Lenton, 2015). Other shifts in environmental stochasticity are associated with sea temperature. For example, the Pacific Decadal Oscillation and North Pacific sea surface temperatures have become more red-shifted in the period from 1900 to 2015 (Huntingford et al., 2013;Boulton and Lenton, 2015). The reddening of such climate variability entails that populations experience prolonged periods of potentially favorable conditions, but most likely also of unfavorable, more extreme conditions (van der Bolt et al., 2018). Prolonged periods of high temperate can skew sex ratios in sea turtles (Spotila et al., 1987), reduce hatchling survival (Spotila and Standora, 1985;Matsuzawa et al., 2002;Glen et al., 2003;Godfrey and Mrosovsky, 2006), and reduce nesting and feeding habitats due to increasing frequency of storms associated with high sea surface temperature (Hawkes et al., 2009;Pike and Stiner, 2007;Pike et al., 2015). More elaborate, species-specific DEB models that include effects of temperature (e.g., Marn et al., 2017aMarn et al., ,b, 2019Stubbs et al., 2020) could form the basis to explore in detail how populations respond to such environmental change. All in all, our findings signify the importance of understanding how life history traits and population responses to environmental change are linked.
Such understanding is a basis for accurate predictions of marine megafauna population responses to environmental perturbations like (over)fishing, and to shifts in the autocorrelation of environmental variables, ultimately contributing toward bending the curve on marine biodiversity loss.

DATA AVAILABILITY STATEMENT
The datasets generated in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://doi.org/10. 6084/m9.figshare.12738749.

AUTHOR CONTRIBUTIONS
IS conceived of the project and took the lead in design, coordination, and writing. All authors brought distinctive expertise to the collaboration and contributed importantly to the ideas represented, as well as through the drafting and revising of the manuscript.

FUNDING
IS acknowledges funding from the Netherlands Organisation for Scientific Research (VIDI grant no. 864.13.005).