Relative contributions of fixed and dynamic heterogeneity to variation in lifetime reproductive success in kestrels ( Falco tinnunculus )

Many species show large variation in lifetime reproductive success (LRS), with a few individuals producing the majority of offspring. This variation can be explained by factors related to individuals (fixed heterogeneity) and stochastic differences in survival and reproduction (dynamic heterogeneity). In this study, we study the relative effects of these processes on the LRS of a Dutch Kestrel population, using three different methods. First, we extended neutral simulations by simulating LRS distributions of populations consisting of groups with increasingly different population parameters. Decomposition of total LRS variance into contributions from fixed and dynamic heterogeneity revealed that the proportion of fixed heterogeneity is probably lower than 10% of the total variance. Secondly, we used sensitivities of the mean and variance in LRS to each parameter to analytically show that it is impossible to get equal contributions of fixed and dynamic heterogeneity when only one parameter differs between groups. Finally, we computed the LRS probability distribution to show that even when all individuals have identical survival and reproduction rates, the variance in LRS is large (females: 27.52, males: 12.99). Although each method has its limitations, they all lead to the conclusion that the majority of the variation in kestrel LRS is caused by dynamic heterogeneity. This large effect of

lifetime reproductive success (LRS). Most individuals have little or no LRS while a few individuals produce the majority of offspring (Clutton-Brock, 1988;Newton, 1989). This variation can be caused by factors related to individuals, for example, genotypic or trait differences that affect survival and reproduction. The LRS variation caused by these deterministic factors is called fixed heterogeneity (e.g., Tuljapurkar, Steiner, & Orzack, 2009). However, stochastic processes also lead to variation in LRS (Caswell, 2009;Newton, 1989;Snyder & Ellner, 2018;Tuljapurkar et al., 2009).
For life histories, a phenotype is simply a set of transition probabilities, consisting of survival probabilities and movements among different reproductive states (e.g., immature, not reproducing, low or high reproductive output) conditional on survival (these are also called switching probabilities; Steiner & Tuljapurkar, 2012). Due to the stochastic nature of these probabilities, individuals with identical life-history phenotypes may have different LRS because they follow different trajectories of reproductive states (between-trajectory variance) and also because individuals that follow the same trajectory of reproductive states die at different ages (withintrajectory variance). Together, the variance in LRS caused by these between-and within-trajectory variances is called dynamic heterogeneity (e.g., Steiner & Tuljapurkar, 2012;Tuljapurkar et al., 2009). Individuals may thus produce a high LRS not because they have some high-quality traits that make them superior, but because they are simply lucky (Snyder & Ellner, 2016. The total variation in LRS is caused by this dynamic and fixed heterogeneity (note that some studies call these patterns individual stochasticity and heterogeneity, respectively, for example, Caswell, 2009Caswell, , 2011van Daalen & Caswell, 2017). Several studies have shown that dynamic heterogeneity constitutes a major part of the total variation in LRS (Caswell, 2011;Jenouvrier, Aubry, Barbraud, Weimerskirch, & Caswell, 2018;Orzack, Steiner, Tuljapurkar, & Thompson, 2011;Snyder & Ellner, 2016Steiner & Tuljapurkar, 2012;Steiner, Tuljapurkar, & Orzack, 2010;Tuljapurkar et al., 2009;Tuljapurkar, Zuo, Coulson, Horvitz, & Gaillard, 2020;van Daalen & Caswell, 2017). In this study, we test whether dynamic heterogeneity is also the dominant driver of LRS variation in a Dutch population of the Common Kestrel (Falco tinnunculus).
In an earlier study on this kestrel population, we tested the effects of several factors which can be expected to influence kestrel LRS (e.g., fledging weight, average laying date of the first egg, age of the mate, quality of the breeding location, LRS of the mother) but only found a significant effect of average laying date and age of the mate (Broekman, 2016). The absence of a significant effect for most factors might be attributed to a large effect of dynamic heterogeneity, which can obscure relations with LRS. Because of this result and the finding in many other studies that dynamic heterogeneity plays a large role in determining LRS variation, we hypothesize that for our kestrel population the majority of LRS variation can also be attributed to dynamic heterogeneity.
To test this hypothesis, we constructed separate population models for males and females, which we used in three different ways to test the importance of dynamic heterogeneity for the LRS distribution in this population. First we simulate LRS distributions of in silico populations of the same size as the observed population and with the same average population parameters, but with increasing differences in demographic values between (groups of) individuals (i.e., a new extension of methods developed by Steiner & Tuljapurkar, 2012;Steiner et al., 2010;Tuljapurkar et al., 2009). We then partition the total variance into variance within groups (dynamic heterogeneity) and variance between groups (fixed heterogeneity). Secondly, we perform sensitivity analysis of the mean and variance in LRS to the population parameters (van Daalen & Caswell, 2017) and use these sensitivities to estimate variance within and between groups of individuals with similar population parameters. Finally, we apply a recently developed method by Tuljapurkar et al. (2020) to compute the probability distribution of producing every possible LRS, as these exact probabilities are more useful descriptors of a highly skewed LRS distributions than the mean and the variance. We then calculate the sensitivities of producing a certain LRS to the population parameters. These methods were performed separately for females and males and for both sexes all methods agree that most of the variation in LRS can be attributed to dynamic heterogeneity, thus confirming our hypothesis.

| Study population and data collection
The Common Kestrel (F. tinnunculus) is a small raptor that is widespread in Europe and is the most common bird of prey in the Netherlands. It is a short-lived species and can start breeding in their second calendar year (Village, 1990). We studied a kestrel population in the Haarlemmermeer (185 km 2 ), which contains Amsterdam Schiphol Airport. This area has on average 56 breeding pairs per year (varying between 36 and 73 breeding pairs per year, Table S1), with clutch sizes between one and eight eggs (Table 1). The LRS in this population ranges between 0 and 42 for females and 0 and 38 for males, with the majority of the offspring produced by a small fraction of the population (4% of the kestrels produce around 30% of the offspring).
Data on this population (Bol, Broekman, Jongejans, & Tuljapurkar, 2020) were collected between 1993 and 2015 and entailed capturing and individually marking of birds, as well as observing the success of all nests (mainly in nest boxes) in the study area. In addition, breeding biology data were gathered, which included, among others, the location, the number of eggs and fledglings and the identity and age of the parents. The identity of the mother was known for 81% of the nests and the identity of the father was known for 29% of the nests (n = 1,194). Furthermore, in a previous study on this population we determined the amount of extra-pair paternities in 109 different nests from 2011 to 2014, including years with both low and high prey abundance and did not find any extra-pair paternities (Broekman, 2016). We can therefore assume that the identified father of a nest is the father off all offspring from that nest.
The main prey item for kestrels in this population are Common Voles (Microtus arvalis) (Village, 1990). Since the reproductive output of the Long-eared Owl (Asio otus) is closely related to vole abundance (Korpimaki, 1992;Sharikov, Volkov, Svhidova, & Buslakov, 2019) and pellet analysis confirmed that Long-eared Owls mainly feed on voles in the Haarlemmermeer (70-80% of all the prey items), we used their reproductive success in the same area and time period (B. J. Bol, unpublished data, 2016, Table S1) to categorize years with low, medium or high prey abundance. Years in which the average number of young was below 1.9 fledglings (1999, 2001, 2003, 2006 and 2013) were categorized as years with a low prey abundance and years in which the average number of young was above 3.0 fledglings (2002, 2004, 2011, 2012 and 2014) were categorized as years with a high prey abundance. All 13 other years were assumed to have an intermediate prey abundance (average of 2.4 fledglings).
Furthermore, to quantify winter harshness, we obtained each year's Hellmann number (KNMI, 2020), which is the sum of all the daily average temperatures below freezing point (without minus sign) in the period November 1 to March 31.

| LRS calculation
We calculated the LRS of kestrels as the sum of the number of fledglings produced in each nest of which they were identified as a parent. This calculation was only performed for kestrels that met the following criteria: (a) they had at least one breeding attempt in the study area, (b) they had no breeding attempt before 1993, (c) they were 2 years old or younger at their first capture, (d) they were not seen outside the study area after having bred in the study area, and (e) they have stopped breeding or died in 2015.
The first criterion was necessary because individuals that did not breed at least once in the study area, could either have died before breeding or have bred outside the study area, so their LRS is unknown. The second criterion involved that all kestrels older than 2 years in 1993 are excluded, because kestrels can start breeding in their second year (Village, 1990) and these kestrels could thus have had a breeding attempt we could not detect. We therefore also excluded individuals that were older than 2 years at their first capture, because they could have had a breeding attempt outside the study area or an undetected breeding attempt inside the study area. Similarly, criterion (d) guaranteed that we did not include individuals for which we could have missed a breeding attempt, because it had a breeding attempt outside the study area. Finally, to meet the last criterion we only included individuals that were not seen for at least 2 consecutive years, following for example, Martinez-Padilla, Vergara, and Fargallo (2017). This makes sure that we did not underestimate the LRS of an individual because it may still be breeding in the future.
This selection yields 301 females and 148 males for which the LRS could be calculated as accurately as possible. However, this selection still includes kestrels that were not observed breeding for one or more years (i.e., are in the unknown reproductive State N, see Population model and parameter estimation). Individuals that were not observed breeding in a year could either be individuals that were not breeding that year or individuals that were breeding that year but were not identified as a parent of a nest. If the latter case is true, we still slightly underestimate the LRS.
For females this is probably not much of a problem, because for most nests we could correctly identify the T A B L E 1 Number of nests that produced a certain number of fledglings, as well as the number of nests for which the mother or father could be identified mother. Therefore, if a female is not identified as a parent in a specific year, there is a high probability that the female was not breeding that year. However, for males this can lead to an underestimation of the real LRS, because for 71% of the nests we could not identify the father. If a male is not identified as a parent in a specific year, there is thus a large probability that it still produced some offspring, which are not incorporated in our LRS calculations. This explains why the average LRS is smaller for males (5.63) than for females (6.92).

| Population model and parameter estimation
To estimate the potential influence of dynamic heterogeneity we constructed separate stage-structured population models for males and females. Based on preliminary model selection (Supporting Information S1), we distinguished five states for both sexes, one immature state containing only first-year kestrels (I) and four adult states: individuals producing 0-2 offspring (A), 3-5 offspring (B), 6-8 offspring (C) or an unknown number of offspring (N, including non-reproductive individuals). For each state we estimated survival (s) and switching probabilities between states (conditional on survival, ψ) using multistate capture-mark-recapture (CMR) analyses (Lebreton & Pradel, 2002). We only incorporated observations of live captures inside the Haarlemmermeer (4,069 observations for females, 3,109 for males). The unknown reproductive output state contains all adult individuals which are not found breeding and therefore includes all individuals which are not captured. The capture probability in all other adult states is thus fixed to 1. Furthermore, kestrels cannot move back from one of the adult states (A, B, C or N) to the immature state (I). Moreover, we assumed reproduction to take place before survival as most kestrels die in the winter following the reproductive period. The main model we studied allowed the survival and switching probabilities to depend on the state. However, we also examined effects of prey abundance, temperature and age. To determine the most parsimonious model for males and females separately, we first performed goodness-of-fit testing on the model with only state effects. We used the parametric bootstrap approach and the median-ĉ approach implemented in program MARK (White & Burnham, 1999), operated using the RMark package (Laake, 2013) in R 3.5.2 (see Supporting Information S1 for more information about the model selection).
Parameter estimates from the most parsimonious model were used to construct a vector with survival estimates (s) and a matrix of switching probabilities (Ψ) conditional on survival, separately for females and males.
For the matrix of switching probabilities ψ IA indicates the probability of switching from State I to State A. These vector and matrix can be combined with the average fertilities (f ) in each state to construct the following population matrix: The fertilities in states A, B and C are calculated separately for females and males as the mean number of offspring produced by all individuals in each state (Table 1). State N consists of both non-reproducing individuals and individuals with an unknown reproductive output, which are all unobserved. The fertility of State N is thus calculated as the mean number of offspring from all non-reproducing adults (=0) and nests in which the mother (for females) or father (for males) was unknown. We assume that nonreproductive adults constitute 5% of the adult kestrel population (B. J. Bol, personal communication).
We use these population models to study the effect of dynamic heterogeneity on LRS variation using three different approaches, which we explain in the next three sections: (a) variance partitioning in simulated LRS distributions, (b) sensitivity analysis and (c) computing the LRS distribution. All these analyses were performed separately for females and males.

| Variance partitioning in simulated LRS distributions
We first quantified the relative influence of fixed and dynamic heterogeneity by simulating populations of the same size and with the same age distribution as the observed population (following e.g., Steiner & Tuljapurkar, 2012;Steiner et al., 2010;Tuljapurkar et al., 2009). The observed population consisted of all kestrels from which the LRS could be calculated.
To obtain simulated populations, we started by simulating 15,000 trajectories of states (I, A, B, C and N) for 15 years (3,000 trajectories for each of the five possible starting states), assuming a first-order Markov Chain with the estimated parameters from matrix Ψ as the transition probabilities and the states as the state variables. When there was an effect of prey abundance on the switching probabilities from matrix Ψ, at each time step we first sampled a matrix of switching probabilities and used the parameters from the sampled matrix to determine the state at the next time step. Out of the 23 years used to estimate switching probabilities, 5 years had a low prey abundance, 13 years had an intermediate prey abundance and 5 years had a high prey abundance. We therefore assumed probabilities of 0.217, 0.565 and 0.217 to sample a matrix of switching probabilities corresponding to years with low, intermediate and high prey abundance, respectively.
To each individual in the observed population we randomly assigned one of these trajectories in the following way (see Box 1, Table 2): we first retained only the trajectories from which the first state corresponds to the state in which the observed individual was first detected.
For each of these trajectories we then calculated the cumulative survival up to the year the observed individual was last detected, using the survival estimates from s and assuming an immature survival of 1 (as we only studied adult kestrels). To accommodate for variation in survival rates due to yearly variation in prey abundance (when the most parsimonious model contained an effect of prey abundance on survival rates), at each time step we determined the survival rate by first sampling a vector of survival probabilities, similar to the sampling of a matrix of switching probabilities. The cumulative survival probabilities were then rescaled to sum to one and used as probabilities to randomly select one trajectory, similar as in for example, Steiner et al. (2010) and Orzack et al. (2011).
This method guarantees that the simulated population has the same age structure as the observed population. This is important because the observed population consist only of a small fraction of all the individuals that were used to estimate the population parameters. This fraction is biased toward older individuals, because individuals are only included in the observed population if they have bred at least once, so the survival rates in the observed population are probably higher than the estimated survival rates. Simulating the LRS distribution directly with the estimated survival and switching probabilities would therefore lead to a younger age structure than the observed population and hence lower LRS in the simulated populations, making the LRS distributions from the observed and simulated population incomparable. Ideally, we would have estimates of survival rates based only on individuals of the observed population. However, as the number of individuals in the observed population is too small to provide reliable estimates, we used the method described above to prevent differences in survival rates between the observed population and the estimated survival rates to affect our results.
After selecting a random trajectory for each individual, we randomly assigned a number of offspring to each state for each individual, using as probability weights the fractions of individuals producing particular numbers of offspring in each state. Because we do not know what number of offspring individuals from the observed population produce when they are in State N and because we wanted the LRS distributions of the observed and simulated populations to be comparable, we assumed that individuals from the simulated populations also produced no offspring when they were in State N (which is different from the population model used to compute the sensitivities and the LRS probability distribution).
Finally, the LRS for each simulated individual is calculated by summing the number of offspring produced each year until the final age of the corresponding kestrel

BOX 1 Numerical example of simulating the LRS of an individual
To illustrate how we simulate the LRS of an individual, we look at the female with ring number 3493066. It was captured for the first time in 1993 as an immature and last observed in 1999, so we assume it survived for 6 years. Ten possible simulated sequences of reproductive states are shown in the first column of Table 2. Survivorship up to age 6 was also calculated for each trajectory and rescaled to sum to one. We use these rescaled survivorships as probability weights to randomly select one trajectory, for example INNNBCB.
For each state in this trajectory we then calculate the number of offspring that are produced. For example, in the fifth year, the individual occurs in State B. Individuals in State B produce three, four or five offspring with probabilities 0.19, 0.34 and 0.48, respectively. Based on these probabilities we assigned four offspring to stage B in the fifth year. Similarly, to the sixth and seventh year, we assigned six and three offspring, respectively. These numbers of offspring are then summed to obtain the LRS, in this case 13. The observed LRS for this female was 18. observed in the field. In this way we derived LRS distributions of 1,000 simulated populations (see also Orzack et al., 2011;Steiner et al., 2010). We tested if the average LRS distribution of these 1,000 simulated populations are different from the observed LRS distribution using a Kolmogorov-Smirnov test.
This method (also called neutral simulations) assumes that every individual has the same phenotype (survival and switching probabilities), so all observed variation is due to dynamic heterogeneity. To assess the potential relative influence of dynamic heterogeneity when variation in LRS is caused by both dynamic and fixed heterogeneity, we extend this method by repeating the simulations with the inclusion of variation in phenotypes between individuals. This involved that before we simulated a trajectory we assigned the simulated individual randomly (with equal probability) to either one of two groups: low or high-quality group. These two groups have population parameters differing from the mean with the same deviation (w), but in opposite directions. Simulations then proceed in the same way as described above.
We studied four different scenarios: (a) fixed heterogeneity in survival, (b) fixed heterogeneity in reproduction, (c) fixed heterogeneity in survival and reproduction with a trade-off between these processes and (d) fixed heterogeneity in survival and reproduction with a positive correlation between them. All these scenarios are simulated with values of w (deviation from the mean) ranging from 0 to 0.15 with steps of 0.01.
To include fixed heterogeneity in survival rates, we perturbed all adult survival rates with the same w. Although survival rates are not used in simulating trajectories, they still change the LRS distributions of simulated populations as they change the cumulative survival probabilities for each trajectory and therefore the chance a trajectory is assigned to a certain individual. To include fixed heterogeneity in reproduction, we perturbed the switching probabilities to the intermediate reproductive state (State B) and the high reproductive state (State C) with value w. We changed the other transition rates in the opposite direction and with the same total change. Differences in switching probabilities change the LRS distributions of simulated populations as these probabilities are used in simulating trajectories.
The total variation in LRS of the simulated populations can be decomposed into variation between the two groups (V between ) and variation between individuals within the same group (V within ) with the following equations: p low is the fraction of individuals within the lowquality group, V low and V high are the variances in the LRS for individuals in the low-and high-quality group, respectively, μ low and μ high are the mean LRS for individuals in the low and high-quality group, respectively, and μ is the mean LRS of the whole population. We also calculate the power of a one-sided t-test to detect a significant difference in the mean LRS between the two groups at the 5% level.

| Sensitivity analysis
In addition to the simulations, we also use sensitivity analyses to study the relative effects of fixed and dynamic heterogeneity. Compared to the simulations, sensitivity analysis offers the advantage that we can approximate for T A B L E 2 Ten simulated trajectories of reproductive states up to Age 6, using the switching probabilities for females ( Note: For each trajectory we calculated the survival up to Age 6 using the survival probabilities for females (Table 3). These survivals were rescaled to sum to 1 to obtain a probability weight for each trajectory. We calculated the number of offspring produced in each year and the resulting LRS of each trajectory by randomly assigning a number of offspring to each State, assuming a 49% probability to produce zero offspring in State A, 19% probability to produce one offspring in State A, 32% probability to produce two offspring in State A, 19% probability to produce three offspring in State B, 34% probability to produce four offspring in State B, 48% probability to produce five offspring in State B, 92% probability to produce six offspring in State C, 7% probability to produce seven offspring in State C and 2% probability to produce eight offspring in State C. We assume no offspring are produced in States I and N. each population parameter how large the difference between the two groups needs to be to get a specific relative contribution of fixed heterogeneity. In addition, calculating sensitivities offers the advantage that it reveals which parameter has the biggest impact on these potential relative influences and less computation is needed than for the simulations. A disadvantage of sensitivity analysis is that it cannot account for yearly variation in population parameters, due to for example, variation in prey abundance. The first step in the sensitivity analysis is therefore to construct a vector of average survival rates and a matrix of average switching probabilities (if the most parsimonious model contained an effect of prey abundance on survival or switching probabilities). To calculate averaged survival and switching probabilities, we weighted the parameters from years with low, intermediate and high prey abundance by the probability that these years occur (i.e., 0.217, 0.565 and 0.217, respectively). These averaged survival vector and matrix with switching probabilities were then used in the equations derived by van Daalen and Caswell (2017) to calculate the sensitivities of the mean and variance of the LRS to the population parameters. Their equations for the mean and variance in LRS are extensions of the equations from Steiner and Tuljapurkar (2012) and incorporate both between-and within-trajectory variance in reproduction. Note that in these equations we assume that individuals can reproduce in the unknown reproductive state (their fertility is given by f N ), which is different from the simulations of the LRS. We can use these sensitivities to estimate the withingroup and between-group variances in the LRS given differences in a specific parameter between the two groups (see Supporting Information S2 for derivations): In which ∂E LRS I ½ ∂s I indicates the sensitivity of the mean LRS to a change in the survival rate of State I. However, this parameter can be substituted for the sensitivity of the mean LRS to any other population parameter to study the effect of a change in other parameters on the change in mean LRS and thus the between-group variance. Equations 4-6 can subsequently be combined to produce an equation that calculates how large a difference in a population parameter between two groups needs to be to get a specific fraction of the total variation in LRS that is between-group variance (i.e., to get a specific relative contribution of fixed heterogeneity): It is important to note that these equations are only valid when the groups have an equal number of individuals and when the relation between the mean or variance in the LRS and the difference in parameter values (w) is linear. The latter assumption is probably only valid for a small w, that is, a small difference in the population parameter between the two groups. We therefore still performed simulations to check whether our calculations with sensitivities lead to the correct conclusions (Supporting Information S3).

| LRS probability distribution
Finally, we computed the probability distribution of producing every possible LRS, assuming the same phenotype for every individual, using the methods described in Tuljapurkar et al. (2020). We started by computing exact probability distributions of reproductive outputs at each state by dividing the number of individuals with a certain reproductive output by the total number of individuals captured in each state (this also includes computing a probability distribution for the unknown reproductive state, which contrasts with the simulations where individuals are assumed to not reproduce in this state). For example, females with four offspring were captured 199 times. In total 592 individuals were captured in State B, so the probability of producing four offspring in State B is 0.34. The mean reproductive outputs in each state, calculated from these probability distributions correspond to the fertilities in the population matrix.
The stage-specific probability distributions can be combined with the survival and switching probability to calculate probability distributions of producing a certain number of offspring at each age. Convolutions of these age-specific probability distributions (weighted by survivorship and mortality) then yield the probability distribution of all possible values of the LRS (see Tuljapurkar et al., 2020 for details). From this probability distribution we can calculate the mean and variance of the LRS, which are similar to the mean and variance calculated with the analytical equations from van Daalen and Caswell (2017). Because this analysis can also not accommodate for yearly variation in population parameters, survival and switching probabilities are averaged over years with low, intermediate and high prey abundance, similar to the sensitivity analysis.
We used the LRS probability distribution to calculate the Gini index, a measure of the inequality of LRS between individuals. A Gini index close to zero indicates similar probabilities of producing each LRS, while a Gini index close to one means that there is a wide variation in the probability of producing a certain LRS and is what we expect for a highly skewed LRS distribution .
Furthermore, we used the LRS probability distribution to calculate the sensitivity of the probability of producing a certain LRS to each population parameter numerically. We added small perturbations (e.g., 10 −9 ) to each population parameter (survival and switching probabilities) and quantified the change in the probability of producing a certain LRS. Division of these changes in probability by the small perturbation yields the sensitivities of the probability to each population parameter. Note that an alternative method of analytically finding sensitivities for the probability of producing a LRS of zero can be devised, with some effort from the methods of Ellner, Childs & Rees (2016, chapters 3,4). However, this method cannot be used for calculating sensitivities for the probabilities of producing a higher LRS. As we prefer to use the same method to calculate sensitivities for the probability of producing every LRS, we calculated all these sensitivities numerically.

| Population model
The most parsimonious population model for females had different survival and switching probabilities for each state and additional effects of years with low prey abundance on the survival rates and of years with high prey abundance on the switching probability to the high reproductive state (State C). The most parsimonious population model for males had effects of state and years with high prey abundance on survival rates and only an effect of state on the switching probabilities. Winter harshness did not significantly affect switching probabilities (the best model including the Hellmann index scored 26.5 ΔQAICc worse than the top model for females and 3.2 for males), while there was some support for winter harshness affecting survival in females, but not males (1.0 ΔQAICc for females and 5.1 for males, see also Supporting Information S1 and Tables S4 and S5). For females the survival rates varied between 0.28 and 0.66, with the highest survival rate in the high reproductive output state and lower survival rates in years with a low prey abundance. For males the survival rates varied between 0.21 and 0.53, with the highest survival rate in the state with the unknown reproductive output and higher survival rates in years with a high prey abundance (Table 3). Furthermore, females had the highest probability of going to the state with the intermediate or unknown reproductive output and a higher probability of going to the state with the high reproductive output in years with high prey abundance. Males had the highest probability to go to the state with intermediate or unknown reproductive output ( Table 4). The average fertilities in states A, B, C and N are 0.83, 4.29, 6.11 and 2.06 for females and 0.93, 4.32, 6.07 and 3.39 for males. These parameters from the most parsimonious models were used to obtain the results described below.

| Simulation results
The mean LRS observed for females was 6.92 with a variance of 31.01 and for males 5.63 with a variance of 28.30, indicating highly skewed distributions. In simulations without fixed heterogeneity, the average mean and variance were 7.08 and 37.13 for females and 5.95 and 21.51 for males. The simulated LRS distributions were not significantly different from the observed distribution for females (KS-test, D = 0.3125, p = .42), while they differed slightly for males (KS-test, D = 0.500, p = .037) (Figure 1).
When fixed heterogeneity increased in the simulations, the means (Figure 2 and S1) and variances ( Figures S2-S3) of the two groups became increasingly different in all scenarios, except the one with a trade-off between survival and reproduction. The population mean LRS only decreased when we assumed that trade-off.
With increasing fixed heterogeneity, the scenarios also showed an increased proportion of the total variance that is between-group variance (Figure 3 and S4), indicating a larger relative effect of fixed heterogeneity (except with the trade-off). This also led to a higher power to detect a difference in LRS between the two groups (Table S6). Nevertheless, the between-group variance never exceeded 10% of the total variance.

| Sensitivity analysis
For both females and males mean LRS is most sensitive to immature survival, followed by survival in the unknown reproductive state. Of the switching probabilities, the probabilities of moving from the immature state and unknown reproductive state are the most important for both sexes, but lower for males (Figures 4 and S5). The variance in LRS is most sensitive to changes in the survival in the immature and unknown reproductive state for both sexes. In addition, for females survival in the intermediate reproductive state (B) is also important.
For the switching probabilities the highest sensitivities are found for the probabilities of moving to the low (A) or high (C) reproductive states, with higher sensitivities for females ( Figures S6-S7).
As we are interested in the relative effects of dynamic and fixed heterogeneity, we implement the sensitivities in Equation 4-6 to estimate the fractions of within-and between-group variance when different population parameters differ between the groups. When we assumed a deviation from the mean (w) of 0.15 (the highest amount of fixed heterogeneity studied in the simulations) and looked at the immature survival rate (the most sensitive population parameter and thus the one that has the greatest ability to increase the between-group variance), the fraction of the total variance that is between-group variance was only 0.036 for females and 0.066 for males (simulations produced similar values, Supporting Information S3). Furthermore, it is unlikely that equality of between-group and within-group variances can be achieved by increasing the difference in only one parameter, because w needs to be larger than 0.5 (i.e., 0.77 for females and 0.57 for males). This would lead to parameter values outside the zero to one range, a result that can be confirmed with simulations (Supporting Information S3).

| LRS probability distribution
Finally, we studied the LRS probability distribution ( Figure S8). The mean and variance in LRS are 2.29 and 27.52 for females and 1.42 and 12.99 for males. (These numbers change to 6.82 and 51.00 for females and 6.36 and 26.66 for males, when we assume that all individuals survive as immatures, as is the case for the observed population which only includes individuals that breed at least once, see Figure 1).
As we saw for the observed and simulated populations, this distribution is highly skewed (even more than for the observed and simulated individuals, because for this probability distributions we also assumed that individuals could die before having a breeding attempt). This highly skewed distribution is indicated by a high Gini-index for both females (0.84) and males (0.87) and a high probability to produce zero offspring (0.73 for females, 0.80 for males). Furthermore, the probability of having a LRS of 10 or higher is only 8% for females and 5% for males.
The probability of producing zero offspring is most sensitive to immature survival for both females (Figure 5) F I G U R E 3 The fraction of the total variance that is attributed to between-group variance (i.e., to fixed heterogeneity) for females. w indicates the deviation of the survival and/or switching probability of the two groups from the mean values. The light-colored bands around the lines indicate the standard deviation in the proportion of total variance that is betweengroup variance. The four panels show different scenarios of fixed heterogeneity [Color figure can be viewed at wileyonlinelibrary.com] F I G U R E 4 The sensitivities of the mean lifetime reproductive success (LRS) to the survival and switching probabilities for females. I denotes the immature state, A denotes the state with low reproductive output (zero to two offspring), B denotes the state with intermediate reproductive output (three to five offspring), C denotes the state with high reproductive output (six to eight offspring), and N denotes the state with unknown reproductive output [Color figure can be viewed at wileyonlinelibrary.com] F I G U R E 5 The sensitivities of the probability of having a lifetime reproductive success (LRS) of zero to the survival and switching probabilities for females. I denotes the immature state, A denotes the state with low reproductive output (zero to two offspring), B denotes the state with intermediate reproductive output (three to five offspring), C denotes the state with high reproductive output (six to eight offspring), and N denotes the state with unknown reproductive output [Color figure can be viewed at wileyonlinelibrary.com] and males ( Figure S9), as increased immature survival will lower the probability of producing zero offspring. Increasing the probability of going from the immature state to the low reproductive state (State A) will slightly increase the probability to produce zero offspring for both sexes, but this effect is low compared to the effect of immature survival. For the probabilities of producing a higher number of offspring, immature survival remains among the most sensitive parameters, but its effect is now positive (Figures S10-S19). Furthermore, switching probabilities from the immature state have the largest effects (compared to switching from other states).

| Evidence for large effects of dynamic heterogeneity on kestrel lifetime reproductive success
Our analyses showed that it is very likely that our hypothesis that the majority of variation in kestrel LRS can be attributed to dynamic heterogeneity is true, as the results from all performed analyses point to this conclusion. We first showed that neutral simulations (without fixed heterogeneity) led to LRS distributions for females similar to the observed LRS distribution, echoing the results of previous studies on other organisms performing neutral simulations (Steiner et al., 2010;Steiner & Tuljapurkar, 2012;Tuljapurkar et al., 2009). Dynamic heterogeneity alone could therefore account for the observed variance in female kestrel LRS, although it is likely that at least part of this variance is due to fixed heterogeneity (Authier, Aubry, & Cam, 2017;Jenouvrier et al., 2018).
Interestingly, variance in LRS of the simulated female LRS populations without fixed heterogeneity (37.13) is even higher than in the observed population (31.01), although some fixed heterogeneity might be expected in the observed population, which could increase the LRS variance. However, we attribute this result to the large variation in the LRS variance between simulations that can arise when the number of individuals in the population are relatively low (301 females). For example, although the average LRS variance of the 1,000 simulations in our study is 37.13, the variance in each simulation ranges from 26.84 to 48.58. Furthermore, a trade-off between survival and reproduction rates can also explain the lower LRS variance in the observed population as fixed heterogeneity in this scenario decreases the total variance ( Figure S2). A similar result has also been found by a previous study (Steiner & Tuljapurkar, 2012), which detected a higher variance in lifespan in simulated populations (without fixed heterogeneity) compared to the observed population.
For male kestrels we found evidence for fixed heterogeneity as there was a slightly significant difference in LRS distribution between the observed and simulated populations (the observed variance in LRS is 28.30, whereas the variance in simulated populations without fixed heterogeneity is only 19.26). However, the simulated LRS distribution for males was still highly skewed, indicating that dynamic heterogeneity leads to large variation in LRS and the difference in results with females might be attributed to the lower fraction of nests from which the identity of the father is known.
Secondly, we showed by simulating LRS distributions of in silico population with increasing differences in populations parameters and decomposing the total variance into contributions from fixed and dynamic heterogeneity, that even with large fixed differences (w = 0.15 for all survival and switching probabilities) the relative effect of fixed heterogeneity is only around 10%.
Thirdly, using sensitivities of the mean and variance in LRS to each population parameters, we could calculate analytically that fixed differences in single population parameters cannot lead to equal contributions of fixed and dynamic heterogeneity and only leads to minor increases in the relative effect of fixed heterogeneity (i.e., approximately a 3.6 or 6.6% increase for females and males, respectively, when we assumed a difference of 0.3 in the most sensitive parameter [immature survival] between two equally sized groups).
Finally, the LRS probability distribution also showed that even when all individuals have the same phenotype, the probability of producing a certain LRS is highly skewed, as is also the case for the observed population. Dynamic heterogeneity can thus produce a large variation in LRS. Because most kestrels die as immatures, immature survival and the probability of moving from the immature state to another state have the largest effect on this distribution. We attribute the differences in mean LRS from the observed and simulated populations compared to the mean LRS calculated from the LRS probability distribution, to the fact that the observed population only includes individuals that had at least one breeding attempt. This is not necessarily the case for the LRS probability distribution, even not when we assume an immature survival of 1 (Figure 1), because kestrels can also die as adults without having had any breeding attempt. Furthermore, differences between the mean LRS from the computed LRS probability distribution (assuming an immature survival of 1) and the observed and simulated mean LRS, arise because for the observed and simulated populations we assume that no offspring are produced when individuals are in the unknown reproductive State N, whereas in the computation of the theoretical LRS distribution, there is still a probability that individuals in State N reproduce.

| Implications of large proportional effects of dynamic heterogeneity for kestrels
The potentially large influence of dynamic heterogeneity on kestrel LRS can have several implications for this population. First of all, as already mentioned in the introduction, large effects of dynamic heterogeneity make it difficult to derive which factors affect LRS as it obscures the relationship between LRS and individual traits (Snyder & Ellner, 2018). In particular, the effects of factors leading to small fixed differences will be difficult to detect. For example, in Broekman (2016) we did not find an effect of fledging weight on kestrel LRS, whereas another study on a Spanish kestrel population found a higher LRS in kestrels with a higher body condition as a fledging (which we assume to be related to fledging weight, Martinez-Padilla et al., 2017). In our population, fledging weight could also be affecting LRS, but it is not detected because it is swamped by the much larger effect of dynamic heterogeneity. Only when a factor leads to large differences in population parameters (see Table S6) or when the sample size is large can the effect of that factor on the LRS be detected (Steiner & Tuljapurkar, 2012).
We also showed that if a factor has opposing effects on survival and reproduction, it will be difficult to detect effects of this factor on LRS. In such cases a variable other than the LRS may be better suited to measure fixed heterogeneity. More generally, our results show that it is important to ask which vital rates are affected by fixed differences.
The large influence of dynamic heterogeneity on kestrel LRS may also imply that most variation in LRS is evolutionary neutral, especially as only a part of the fixed heterogeneity will be due to genetic differences (Hartemink & Caswell, 2018;van Daalen & Caswell, 2017). The consequence would be that a high selection pressure on traits is absent and it can contribute to a low heritability of traits closely related to fitness (Steiner & Tuljapurkar, 2012), which is a common observation in natural populations of many different species (e.g., Kruuk et al., 2000;McCleery et al., 2004;Merila & Sheldon, 2000).
However, Snyder and Ellner (2018) showed that whether a high selection pressure on traits is indeed absent depends on the population size, as in a large population variation due to dynamic heterogeneity averages out over individuals with similar traits. We estimate that our population contains around 120 mature individuals, based on the average yearly number of breeding pairs and a 5% proportion of non-breeding adults. Calculating whether this population size is large enough to limit the effects of dynamic heterogeneity on selection requires estimates of a continuously varying trait and its effect on the expected yearly offspring number. We could then apply the equations from Snyder and Ellner (2018). However, performing these analyses is beyond the scope of our study. Furthermore, we did not look at specific trait values on which selection could act. We can therefore not conclude that a high selection pressure is absent, but the large relative influence of dynamic heterogeneity makes this at least a possibility.

| Dynamic vs fixed heterogeneity in the literature
Whether dynamic or fixed heterogeneity plays a dominant role in determining the total variance is debated by many studies. Some studies demonstrate the importance of fixed heterogeneity (e.g., Aubry, Cam, Koons, Monnat, & Pavard, 2011;Cam et al., 2013;Chambert, Rotella, Higgs, & Garrott, 2013;Oosthuizen et al., 2019;Paterson, Rotella, Link, & Garrott, 2018), while others highlight the importance of dynamic heterogeneity (e.g., Caswell, 2011;Jenouvrier et al., 2018;Snyder & Ellner, 2018;Steiner & Tuljapurkar, 2012;Tuljapurkar et al., 2009Tuljapurkar et al., , 2020van Daalen & Caswell, 2017). Our study complements the findings of a few recent studies that partition the variance in LRS or longevity into contributions from dynamic and fixed heterogeneity (using different methods) and found dynamic heterogeneity to be much more important than fixed heterogeneity (Hartemink & Caswell, 2018;Hartemink, Missov, & Caswell, 2017;Jenouvrier et al., 2018;Snyder & Ellner, 2018). Snyder and Ellner (2018) derived analytical equations to partition the total variance in LRS into contributions from trait variation and luck (stochastic variation) and found that luck usually dominates the variance of LRS. These analytical equations not only allow to assess the relative effect of dynamic heterogeneity, but also to predict when we expect mainly dynamic or fixed heterogeneity to be important. Jenouvrier et al. (2018) classified individuals of the Southern Fulmar (Fulmarus glacialoides) in different groups based on finite mixture models that assign individuals randomly to a hidden state. They estimated vital rates of each state for each group which was used to calculate the variance in LRS and longevity within the groups of individuals (dynamic heterogeneity) and between the groups of individuals (fixed heterogeneity). Seventy-eight percent of the total variance in LRS appeared to be due to dynamic heterogeneity. In addition, Hartemink et al. (2017) and Hartemink and Caswell (2018) used a comparable approach to study variation in longevity in humans and nine different species of nematodes and insects. They also found small effects of fixed heterogeneity. We expect that applying these analytical approaches to our kestrel population would lead to the same conclusions as our simulations and sensitivity analysis.

| Methodological comparisons
In our analyses, we used new extensions of methods from Tuljapurkar et al. (2009), Steiner andTuljapurkar (2012) and van Daalen and Caswell (2017) to assess the effects of dynamic and fixed heterogeneity. Other studies have used different methods to reach similar results for other populations. For example, analytical solutions to partition the total variance in LRS into contributions from fixed and dynamic heterogeneity are provided by for example, Hartemink et al. (2017), Jenouvrier et al. (2018, Hartemink and Caswell (2018) and Snyder and Ellner (2018). However, we preferred the use of simulations as they are intuitive, stay close to the empirical data and are easy to implement. In our analyses we split the population into two groups, but we can easily extend this to multiple groups (see Supporting Information S4). It is therefore also possible to combine our approach with for example, the one from Jenouvrier et al. (2018), who first split individuals into multiple groups based on a finite mixture model and then run the simulations and decompose the variance as described here. To check the robustness of our findings we also applied other methods, that is, analytical calculation with sensitivities and computing the LRS probability distribution. The results of these complementary analyses were in line with the findings of our simulation results.
Another approach to assess if dynamic or fixed heterogeneity are important is to construct models with fixed heterogeneity as a random individual intercept and use multi-model interference to select the most parsimonious model (e.g., Aubry et al., 2011;Aubry, Koons, Monnat, & Cam, 2009;Cam et al., 2013;Chambert et al., 2013;Gimenez, Cam, & Gaillard, 2018;Paterson et al., 2018). However, this approach has several disadvantages which prevented us from using it. First of all, model selection is based on a trade-off between bias and precision in the parameters of the assumed model, which are determined by comparing model results to the same data to which the model was fitted. No different data are therefore used to validate the model. A better way to study the influence of a fixed factor is to split the dataset in test and training datasets, estimate model parameters based on the training dataset and look at the deviance of this model output with the test dataset. Secondly, individual random effects are difficult to interpret biologically. Finally, random effects may reflect factors not in the model and may be sensitive to outliers (Knape, Jonzen, Skold, Kikkawa, & McCallum, 2011). For example, Plard et al. (2015) showed that after accounting for mass and age structure in a population of roe deer, random intercept effects only had a very small effect.

| Limitations of our methods
As all our methods point to a large effect of dynamic heterogeneity on kestrel populations, we can reliably conclude that most of the variation in kestrel LRS is due to dynamic heterogeneity. However, we want to note a few limitations of our analyses. First of all, we simulated LRS distributions with fixed heterogeneity by adding perturbations to the population parameters, so these fixed differences were not observed. If data on real fixed differences are available (e.g., genetic differences), it is better to use these data to group individuals and estimate parameters of these groups. This leads to estimation of real relative effects of dynamic and fixed heterogeneity, instead of potential effects. Secondly, as already mentioned before, our analytical calculations of the proportions of fixed heterogeneity with the sensitivities, assume a linear relationship between the differences in population parameters and the mean and variance in LRS. Nevertheless, simulations showed that this assumption does not severely bias our result. Thirdly, the sensitivity calculations assume an equal number of individuals in both the low-and high-quality group, which might not be realistic. However, the equations can easily be changed to account for other proportions of individuals in both groups (see Equations S4 and S7 in Supporting Information S2). Finally, we lack data to calculate the consequences of large dynamic heterogeneity on trait selection. Although it was not our goal to estimate the effect of dynamic heterogeneity on trait selection, this could produce interesting results, as the lack of a high selection pressure could have negative consequences for the persistence of populations today as climate and environmental changes may demand rapid adaptation (Urban, 2015;Vitousek, Mooney, Lubchenco, & Melillo, 1997).

| Future directions
Our results demonstrate large effects of dynamic heterogeneity in a Dutch kestrel population. We expect similar findings in many other species, as a high entropy of reproductive trajectories and small persistence in states, similar to our population, are common (see Supporting Information S5 and Tuljapurkar et al., 2009). Even in species with large fixed differences in reproduction and hence a low entropy and high persistence in states, stochasticity in survival can still lead to large effects of dynamic heterogeneity. For example, van Daalen and Caswell (2015) found that variance in human LRS within trajectories are larger than between trajectory. Using our analyses, or methods from for example, Snyder and Ellner (2018), it is now possible to test the potential influence of dynamic heterogeneity in many animals, using for example the COMADRE animal matrix database (Salguero-Gomez et al., 2016).
Simulations with fixed differences can also be used to study other scenarios of fixed differences between groups and with more groups (see Supporting Information S4). However, the most interesting next step would be to repeat the analyses with individuals assigned to different groups based on real differences that can be related to fixed heterogeneity such as genetic differences. This would yield indisputable results about the relative effects of dynamic and fixed heterogeneity. White, G. C., & Burnham, K. P. (1999)