A functional trait approach to identifying life history patterns in stochastic environments

Abstract Temporal variation in demographic processes can greatly impact population dynamics. Perturbations of statistical coefficients that describe demographic rates within matrix models have, for example, revealed that stochastic population growth rates (log(λ s)) of fast life histories are more sensitive to temporal autocorrelation of environmental conditions than those of slow life histories. Yet, we know little about the mechanisms that drive such patterns. Here, we used a mechanistic, functional trait approach to examine the functional pathways by which a typical fast life history species, the macrodetrivore Orchestia gammarellus, and a typical slow life history species, the reef manta ray Manta alfredi, differ in their sensitivity to environmental autocorrelation if (a) growth and reproduction are described mechanistically by functional traits that adhere to the principle of energy conservation, and if (b) demographic variation is determined by temporal autocorrelation in food conditions. Opposite to previous findings, we found that O. gammarellus log(λ s) was most sensitive to the frequency of good food conditions, likely because reproduction traits, which directly impact population growth, were most influential to log(λ s). Manta alfredi log(λs) was instead most sensitive to temporal autocorrelation, likely because growth parameters, which impact population growth indirectly, were most influential to log(λ s). This differential sensitivity to functional traits likely also explains why we found that O. gammarellus mean body size decreased (due to increased reproduction) but M. alfredi mean body size increased (due to increased individual growth) as food conditions became more favorable. Increasing demographic stochasticity under constant food conditions decreased O. gammarellus mean body size and increased log(λ s) due to increased reproduction, whereas M. alfredi mean body and log(λ s) decreased, likely due to decreased individual growth. Our findings signify the importance of integrating functional traits into demographic models as this provides mechanistic understanding of how environmental and demographic stochasticity affects population dynamics in stochastic environments.

Interactions between environmental conditions and demography arise because variation in environmental conditions translates into variation in the demographic rates of survival, growth and reproduction, thereby affecting population sizes (Boyce, Haridas, & Lee, 2006;Coulson et al., 2011). The dynamics of such change have been studied extensively in stochastic population analyses where temporal variation in environmental conditions lacks temporal autocorrelation (Getz & Haight, 1989;Lande, Saether, & Engen, 1997).
Results from these analyses reveal that long-lived species with slow life histories, characterized by high survival and low fecundity (Gaillard et al., 1989), are generally buffered from increased environmental variation (Morris et al., 2008, Morris et al., 2011Dalgleish, Koons, & Adler, 2010;Saether et al., 2013; but see Jongejans, Kroon, Tuljapurkar, & Shea, 2010;McDonald et al., 2017). Short-lived species with fast life histories characterized by high fecundity but low survival (Gaillard et al., 1989), on the other hand, are more likely to show increasing fluctuations in population sizes with increasing environmental variation (McDonald et al., 2017;Morris et al., 2008).
In nature, environmental fluctuations usually show positive autocorrelation over time (Ariño & Pimm, 1995;Halley, 1996;Inchausti & Halley, 2002). Yet, we are only beginning to unravel whether links between life history characteristics and population responses to environmental fluctuations identified for uncorrelated environments (Dalgleish et al., 2010;Morris et al., 2008) hold in biologically more realistic, autocorrelated environments. For example, in mite microcosm experiments, life stage responses to environmental perturbations (harvesting) differ between autocorrelated and uncorrelated environments (Cameron et al., 2016;Smallegange & Ens, 2018). Theoretical studies have shown that a change in the serial correlation of demographic rates may increase or decrease population growth rate depending on the structure of the life history (Tuljapurkar, Gaillard, & Coulson, 2009). Using a cross-taxonomical approach, Paniw et al. (2018) most recently illustrated that species characterized by fast life histories exhibit the highest sensitivities to simulated autocorrelation in demographic rates, whereas slow life history species were less sensitive to temporal autocorrelation.
With predicted global changes in environmental patterning (García-Carreras & Reuman, 2011), it is urgent to gain a mechanistic understanding of the life history processes that result in these distinct demographic responses between slow and fast life histories to temporal autocorrelation in environmental conditions. The demographic approaches used to date to map out life history patterns across stochastic environments, however, lack a mechanistic representation of the biological processes that give rise to observed demographic variation (Salguero-Gómez, Violle, Gimenez, & Childs, 2018), as the state variables used in the models (e.g., body size, developmental stage, age) are too phenomenological to explore underlying mechanisms of variation (Salguero-Gómez, 2017).
Mechanistic trait-based models can provide functional insights into the links between individual life history, the environment and population dynamics. Recently, demographic functions describing growth and reproduction have been derived from a simple dynamic energy budget (DEB) growth model, and incorporated into an integral projection model (IPM) (Easterling, Ellner, & Dixon, 2000) to describe survival, growth and reproduction in relation to the feeding environment (Smallegange, Caswell, Toorians, & Roos, 2017).
Here we used such DEB-IPMs to examine the functional pathways by which a typical fast and a typical slow life history species differ in their response to shifts in temporal autocorrelation when their individual growth and reproduction are described by functional, life history traits. We conducted a field study on a barrier island salt marsh to parameterize a DEB-IPM for the fast-reproducing and short-lived macrodetritivore talitrid beach hopper (Orchestia gammarellus) and used a published DEB-IPM of the slow-reproducing, long-lived reef manta ray (Manta alfredi) (Smallegange et al., 2017).
This study is also a demonstration of how a unified method like the DEB-IPM can be applied to species with extremely different life cycles. For each DEB-IPM, we carried out stochastic simulations to assess if the fast (or slow) life history species showed high (or low) sensitivity of the stochastic population growth rate (log(λ s )) to shifts in temporal autocorrelation in food availability (Morris et al., 2008;Paniw et al., 2018). Because environmental effects can affect individual growth and thereby the size of individuals (DeLong et al., 2015;Woodward et al., 2005), we also explored patterns in mean body size across the variable environments. In the stochastic demographic model, the temporal sequence of good and bad food environments was driven by a Markovian process that governs the serial correlation of environment states. In order to cover a wide range of stochastic environments, we varied this serial correlation from blue to white to red noise color (corresponding to a negative first-order autocorrelation, no autocorrelation, and a positive first-order autocorrelation of the temporal sequence, respectively), but also varied the frequency at which the good environment occurred from zero to one. We also conducted a perturbation analysis for each species to assess which life history trait had the strongest effect on the stochastic population growth rate, as this mediates population-level responses to temporal autocorrelation (Franco & Silvertown, 1996Tuljapurkar & Haridas, 2006). Finally, we assessed to what extent patterns in log(λ s ), mean body size and elasticities depended on demographic stochasticity (inherent randomness in demographic processes (Lande et al., 2003), for example, when individuals of the same size and in the same environment, independently of each other, have good or bad luck in their feeding experience), as demographic stochasticity, population structure and temporal autocorrelation can interact to shape population dynamics (Vindenes & Engen, 2017).

| Brief model description
If a female survives from time t to time t + 1, she grows from length L to length L′ following a von Bertalanffy growth curve. If a surviving female is an adult, she also produces eggs. These events are captured by four fundamental functions: the survival function, S(L(t)), the growth function, G(L′,L(t)), the reproduction function, R(L(t)), and the parent-offspring function, D(L′,L(t)), which describes the association between the body length of the parent L and offspring length L′. Denoting the number of females at time t by N(L,t) means that the dynamics of this body length number distribution from time t to t + 1 can be written as: where the closed interval Ω denotes the length domain. The survival function S(L(t)) in Equation 1 is the probability that an individual of length L survives from time t to t + 1: where E(Y) is the expected feeding level, which ranges from zero (empty gut) to one (full gut). Here, individuals that experience E(Y) = 1 can be assumed to always have a full gut. The parameter L m is the maximum body length under conditions of unlimited resource, E(Y) = 1, and κ is the fraction of ingested energy allocated to respiration. 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) < κ).
The function G(L′,L(t)) is the probability that an individual of body length L at time t grows to length L' at t + 1, conditional on survival, and, following common practice (Easterling et al., 2000), follows a Gaussian distribution: with and where ⋅ r B is known as the von Bertalanffy growth rate and σ(Y) is the standard deviation of the expected feeding level. Note that we assumed that individuals can shrink under starvation conditions. The reproduction function R(L(t)) gives the number of offspring where E L b L (t) 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.

| Field study to collect O. gammarellus life history data
The litter feeder O. gammarellus is a terrestrial amphipod belonging to the family Talitridae ( Figure 1a). It is a bioturbating soil animal that typically lives in the litter and soil surface layer of tall, salt marsh intermediate-and late-successional-stage vegetation (Schrama, Berg, & Olff, 2012). At our field site, the salt marsh of the barrier Island of Schiermonnikoog, the Netherlands (53°28′N and 6°12′E), the animals were found between 1.30 and 1.54 m above mean sea level.
They tolerate regular submersion in saltwater due to their osmoregulatory capacity (Morrit, 1989). As is typical of most soil fauna, they are sensitive to high temperatures and drought (Moore & Francis, 1986a). In temperate regions, this species has a yearly life cycle with two reproductive periods, the first at mid-June/early July and the second in early-mid-August ( Figure 1a). Both cohorts overwinter in the subadult or adult stage and reproduce in the summer the year otherwise , after they were born. Mortality mainly occurs after the second reproductive peak. Talitridae do not lay eggs in the soil but carry eggs and embryos in a brood pouch between their legs, in which an average female carries about 17 eggs. Populations are typically female biased (Moore & Francis, 1986b), with on average four times as many females as males (M. P. Berg, unpublished data).  (Wildish, 1969). Juveniles are born with four podomeres and females carrying eggs typically have 12 or more podomeres, while the oldest individuals have 21 podomeres as a maximum.

| DEB-IPM parameterization for O. gammarellus
We used O. gammarellus life history data collected from our field site (see above) to calculate input values for our model analyses. Mean body length at birth was measured at 3.80 ± 0.37 SD (standard deviation) mm. Females matured at a mean body length of F I G U R E 1 (a) Schematic representation of the life cycle of Orchestia gammarellus. The first cohort of juveniles is born between mid-June and early July (arrow 1), depending on the temperature conditions in late winter and spring. Probably, the same females produce a second cohort of juveniles about 6 weeks later in early August (arrow 2). Then, the females that have reproduced die in October (cross) 7.28 mm ± 0.50 SD mm and can reach a maximum length of 15.61 mm.
Based on these observations, we set L b = 3.80 mm, L p = 7.28 mm and L m = 15.61 mm. There are two reproduction peaks, one in mid-June/ early July and the second one in early-mid-August; the highest number of offspring produced by a female in a single clutch during one of these peaks was 32. In the laboratory, under highly favorable conditions, the time between mating and releasing young can be as short as 4 weeks (Morritt & Stevenson, 1993 However, our perturbation analysis (see Section 3.3, Figure 4) revealed that L p and L m were most influential to the log stochastic population growth rate, from which we infer that model outcomes are highly insensitive to the specific values of the von Bertalanffy growth rate and mortality rate. Finally, we assumed

| DEB-IPM parameterization for the reef manta ray
The reef manta ray M. alfredi (Mobulidae) (Figure 1b) is one of the largest rays in the world and is distributed worldwide in tropical and subtropical waters. They are nonmigratory and live close to coasts, reefs or islands, where they aggregate at cleaning stations and areas where they feed on zooplankton (Marshall et al., 2011). We chose to use the DEB-IPM with the additional assumption that individuals cannot shrink under starvation conditions (as shrinking is less likely to occur in vertebrates than in invertebrates) (Smallegange et al., 2017). We use published life history data on female reef manta rays obtained from a stable population off the coasts of Yaeyama Islands, Japan, which population growth rate is estimated at 1.02-1.03 per year (Kashiwagi, 2014). We use the disk width, that is, the distance between the two pectoral fin tips, as the measure for body length. Length at birth is measured at 130 cm and females mature at about 10 years of age at a minimum length of 380 cm (Kashiwagi, 2014) and live at least 40 years, reaching a maximum length of 550 cm (Marshall et al., 2011). On average, adult females produce one pup every 2 years (Kashiwagi, 2014), but this can be as high as one pup every year (Marshall et al., 2011). Based on these observations, we set L b = 130 cm, L p = 380 cm, L m = 550 cm and R m = 1/year. The survival rate of juveniles and adults is estimated at 0.95 (Kashiwagi, 2014) and we set the mortality rate constant at = − log 0.95 = 0.05∕year. The von Bertalanffy growth rate is estimated for females at ⋅ r B = 0.18 (Kashiwagi, 2014). We assumed For the starvation condition, we assumed the commonly used value of κ = 0.80 (Add-my-pet, 2016). All data are summarized in Table 1; note that the time step t of the DEB-IPM is 1 year for the reef manta ray.

| 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. The temporal frequency at which the good environment occurs 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, above and below the feeding level where we predict population equilibrium at λ = 1, we defined E Y low as the expected feeding level as- ) was generated (Smallegange et al., 2014). This sequence determines the environment state, and hence the feeding level

| Perturbation analysis
We conducted a perturbation analysis to examine the elasticity of λ s to perturbation of each of the model parameters L b , L p , L m , R m , ⋅ r B , and µ (Table 1) to identify which life history parameters under which stochastic regimes were most influential to the long-run, log stochastic population growth rate log(λ s ). For this, we perturbed each parameter by 1% and calculated the elasticity of log(λ s ) to each model parameter. Again, we varied the standard deviation of the expected feeding level, σ(Y), and ran the perturbation analysis first with σ(Y) = 0.1, and then with σ(Y) = 0.5.
We excluded the parameter κ, because, apart from occurring in the starvation condition (Equation 2), κ is also related to L m (which is mathematically proportional to κ) and R m (which is mathematically proportional to [1 − κ]) (Kooijman & Metz, 1984). A perturbation analysis of κ is, therefore, mathematically challenging and biological interpretation of its results difficult as effects of variation in κ will carry over to affect survival, growth as well as reproduction. Finally, κ is the only parameter that cannot be measured directly on individuals, which would make it impractical as an indicator of population performance.

| D ISCUSS I ON
Our goal was to examine the functional pathways by which a typical fast and a typical slow life history species differ in their response to shifts in temporal autocorrelation when their individual growth and reproduction are described by functional, life history traits. We chose our study species such that this study is also a demonstration of how a unified method like the DEB-IPM can be applied to extremely different species. The reef manta ray M. alfredi (slow life history) consistently responded to a shift from red to blue noise with a decrease in stochastic population growth rate, regardless of the frequency at which good environments occurred. In contrast, the stochastic population growth rate of O. gammarellus (fast life history) either decreased or increased in response to a shift from red to blue noise, depending on whether the good environment frequency was high or low. More significantly, O. gammarellus stochastic population growth rate was more sensitive to shifts in good environment frequency than to shifts in temporal autocorrelation of food conditions: the higher the frequency value, the higher was the stochastic population growth rate. Qualitatively equivalent patterns emerged when the same stochastic analysis was applied to another fast life history species, the bulb mite Rhizoglyphus robini (Smallegange & Ens, 2018).

| Life history patterns across stochastic environments
Our findings on life history speed sensitivity to shifts in environmental autocorrelation from a mechanistic analysis are opposite to F I G U R E 4 Perturbation analysis results when variability in expected feeding level is low, (a) and (b), or high, (c) and (d). In case of Orchestia gammarellus, the stochastic population growth rate (log(λ s )) was nearly always most sensitive to perturbation of length at puberty (L p ) (light gray areas) (a) and (c), except for a small region at very low values of q, the probability of moving from the bad to the good environment, and low variability in expected feeding level, where log(λ s ) was most sensitive to perturbation of maximum length (L m ) (dark gray areas) (a). In case of Manta alfredi, log(λ s ), was either most sensitive to perturbation of L m , length at birth (L b ) or maximum reproduction rate (R m ) (b) and ( the recent outcome of the phenomenological, cross-taxonomical analysis by Paniw et al. (2018). Although we only have stochastic functional trait analysis outcomes for three species (this study; Smallegange & Ens, 2018), we find there is merit in discussing potential reasons for why the two sets of results differ. Such a discussion is also urgent in light of the persistent problems that exist in the construction of standard matrix population models (Kendall et al., 2019), such that Kendall et al. (2019) advocate to be cautious when interpreting results from, for example, cross-taxonomical analyses that use such models. To start, it is important to note that the finding by Paniw et al. (2018) that fast life histories were on average more sensitive to environmental autocorrelation than slow species only explained 50% of the variance in sensitivity to autocorrelation.
Our study species could therefore simply fall into the category of species that drive the unexplained variance. However, they might not in which case other explanations should be considered that focus on the difference in model approach. Standard stochastic analyses of matrix models perturb demographic rates directly to explore population consequences (Boyce et al., 2006;Morris et al., 2008;Paniw et al., 2018). Because demographic rates in most such matrix models are statistical descriptors (with a mean and variance) of empirical observations (but see Klanjscek, Caswell, Neubert, & Nisbet, 2006), they lack a mechanistic underpinning of how environmental variation impacts demography through individual life history trajectories. An automatic assumption of such an model approach is that statistical descriptors of demographic rates apply to all explored environments (elsewhere we explain what problems can arise if this assumption is not justified; Smallegange et al., 2017).
In a mechanistic approach, like ours, expected feeding level within a stochastic food time series defines the state of the current environment, which subsequently determines an individual's growth and reproduction as dependent on its body size, and an energy allocation trade-off (Smallegange et al., 2017). This presence versus absence of a mechanistic description of demography could explain the contrasting outcomes of our stochastic functional trait analysis and of standard stochastic analyses (Paniw et al., 2018). Further analyses are clearly warranted to confirm or reject this hypothesis, but we may need to be cautious with regard to the assumed mechanisms of how environmental variation affects demography. It also signifies the potential importance of integrating functional traits into demographic models . The fact that in our functional trait approach species only differed in life history trait values, and not in model structure as is often the case in standard stochastic analyses (Paniw et al., 2018), should facilitate future cross-taxonomical comparative studies to evaluate which functional traits best describe and predict life history patterns and population dynamics across variable environments.
Interpreting a population's response to shifts in environmental patterning requires detailed knowledge of the variation in underlying demographic rates (Ezard & Coulson, 2010;Morris & Doak, 2004

| Body size as population performance indicator
Understanding how populations and communities respond to shifts in environmental patterning and how they can transition from stability to decline is an ongoing challenge in ecology (Carpenter et al., 2011;Folke et al., 2004;Scheffer et al., 2009 (Schrama, Boheemen, Olff, & Berg, 2015).
It will be interesting to investigate how large a plastic or genetic shift in life history trait values should be to counter any negative effects of environmental change signaled by a shifts in body size distribution.

| Demographic stochasticity
Finally, we investigated the role of demographic stochasticity in all responses to shifts in environmental quality patterning.
Demographic stochasticity has for a long time been ignored as its impact on the long-term population growth rate levels off with population size (May, 1973). Recent insights, however, highlight how demographic variance, population structure and temporal autocorrelation of environmental conditions can interact to shape population dynamics (Vindenes & Engen, 2017 Roos & Persson, 2013), the resilience of food webs to disturbance (Woodward et al., 2005), and even the strength of trophic cascades (DeLong et al., 2015), investigating the life history mechanisms behind how species of different life history respond to shifts in demographic stochasticity plays will be a valuable effort.

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

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R S' CO NTR I B UTI O N
Isabel Smallegange: conception and design, construction of the population model and performance of the simulations, interpretation of the results, writing of the paper. Matty Berg: conception and design, collection and interpretation of field data, contributed to the interpretation of the results, the writing of the paper, review, and edition of the manuscript prior to submission.