Optimal reproductive phenology under size‐dependent cannibalism

Abstract Intra‐cohort cannibalism is an example of a size‐mediated priority effect. If early life stages cannibalize slightly smaller individuals, then parents face a trade‐off between breeding at the best time for larval growth or development and predation risk from offspring born earlier. This game‐theoretic situation among parents may drive adaptive reproductive phenology toward earlier breeding. However, it is not straightforward to quantify how cannibalism affects seasonal egg fitness or to distinguish emergent breeding phenology from alternative adaptive drivers. Here, we devise an age‐structured game‐theoretic mathematical model to find evolutionary stable breeding phenologies. We predict how size‐dependent cannibalism acting on eggs, larvae, or both changes emergent breeding phenology and find that breeding under inter‐cohort cannibalism occurs earlier than the optimal match to environmental conditions. We show that emergent breeding phenology patterns at the level of the population are sensitive to the ontogeny of cannibalism, that is, which life stage is subject to cannibalism. This suggests that the nature of cannibalism among early life stages is a potential driver of the diversity of reproductive phenologies seen across taxa and may be a contributing factor in situations where breeding occurs earlier than expected from environmental conditions.

& Urban, 2016; Rasmussen, Van Allen, & Rudolf, 2014;Salamolard, Butet, Leroux, & Bretagnolle, 2000;Shulman et al., 1983;Sniegula, Golab, & Johansson, 2019). For instance, a large difference in arrival times of nymphs of two dragonfly species causes the exclusion of a late arrival species (Rasmussen et al., 2014). This is known as priority effects, which emphasize how the sequence of breeding or other phenologies determine interaction strength. Priority effects are often size-mediated (Rasmussen et al., 2014;Sniegula et al., 2019), as organisms emerge small and grow larger while their role as predators and prey shifts (Nosaka, Katayama, & Kishida, 2015). An early start can be a large advantage both in competitive and in predatorprey interactions, but is traded against the match with key resources during ontogeny (Durant, Hjermann, Ottersen, & Stenseth, 2007).
If there is a seasonal peak in the environmental growth or survival conditions for eggs and larvae, cannibalism during early life stages introduces a trade-off for parents between breeding at a favorable time in the season and the added predation risk caused by larger individuals born earlier. Similarly, those parents breeding earlier in the seasonal cycle may benefit by increased survival in their progeny, since they can feed on the younger and smaller offspring arriving later in the season, while availability to other resources in the environment select against earlier breeding. This mechanism may induce seemingly maladaptive breeding phenologies, because the cannibalistic mortality is less evident to us than the environmental drivers. As an example, a recent study pointed out that the Atlantic bluefin tuna spawn surprisingly early in the Mediterranean Sea (Reglero et al., 2018) (Figure 1). The tuna spawn at times when temperatures are lower than the optimum for growth and development for eggs and larvae. This may of course be due to other constraints, such as a food limitation, but the highly cannibalistic nature of the larvae of these tuna (Reglero, Urtizberea, Torres, Alemany, & Fiksen, 2011) may also contribute to this phenology. Other examples of observations where cannibalism may be shaping breeding phenology are the spatiotemporal split of larval load in the fire salamander (Segev et al., 2011), or the early hatching of damselfly (De Block et al., 2005). Such observations suggest that cannibalism can play a major role for the optimal phenology of a wide variety of species. The frequency-dependent nature of the problem requires a noncooperative game-theoretic approach to find the evolutionary stable strategy (ESS) (Iwasa, Odendaal, Murphy, Ehrlich, & Launer, 1983;Krivan, Cressman, & Schneider, 2008) under size-dependent cannibalism.
Here, we devise a theoretical approach to investigate how size-dependent cannibalism can shift breeding phenology in a seasonal environmental cycle. Our focus is on size-mediated priority effects on reproductive phenology, such as the effects of (a) intensity of cannibalism, (b) the ontogeny of cannibalism, and (c) duration of breeding season. We combine age-structured population dynamics with game-theoretic decision making to predict an evolutionary stable breeding phenology. Our goal is to achieve general predictions about the role of intensity and ontogeny of cannibalism, and the duration of the breeding season for breeding phenology. To keep the analysis simple, we focus only on the mortality part of cannibalism as in Diekmann et al. (1986); Hastings and Costantino (1987); Hastings (1987); and Olsson and Andersen (2018). This means that we do not consider resource limitation in the growth of larvae explicitly, and there is no competition for food resources involved in the model, only the extra death risk involved for later born offspring.
In relation to the bluefin tuna example, this is equivalent to assuming temperature is the main environmental driver for growth and development and that larvae will be able to find enough food even without intra-cohort cannibalism, but that they will consume smaller F I G U R E 1 A real-world example of potential cannibalistic drive toward early spawning in the highly cannibalistic Atlantic bluefin tuna larvae in the Mediterranean Sea. This shows the observed spawning phenology (larvae found in the field, bars) and the theoretical fitness of an egg (blue dashed line, scaled to the maximum value) calculated by temperature-dependent egg development time, larvae growth rate, and mortality born at any given day of the year. The red line shows the seasonal temperature cycle, which is an important driver egg and larval development rate and fitness. The figure is modified from Reglero et al. (2018) con-specifics at encounter. This gives a conservative assessment of the benefits of early breeding from a cannibalistic drive.
A common starting point to understand the timing of breeding or spawning in seasonal environments is the "match-mismatch" theory, where parents are selected to match their egg laying to the best environmental conditions for their offspring (Cushing, 1990;Durant et al., 2007). In theoretical models, this is often inferred as a seasonal peak in food availability, which is also the best match for the critical period of the young. We define the seasonal environment from some abiotic variables (e.g., temperature) that determine hatching success at time t. The optimal breeding time is to match the best environment for hatching if there is no cannibalism (Figure 2). Then, we add cannibalism among early life stages and game-theoretic decision making of breeding time among parents. If cannibalism is sufficiently strong, it can lead to deviations from the expectation of all offspring appearing at the optimal environmental conditions for hatching. The intensity of cannibalism is quantified by emerging survivorship through the egg and larval stages, depending on the number of larger individuals in the cohort. Here, we discuss three general size-dependent, ontogenetic cannibalistic modes: (a) Egg are cannibalised by larvae (e.g., anuran, Crump, 1983; and Tribolium (flour beetle), Hastings, 1987); (b) larvae are cannibalised by larger larvae (e.g., damselfly, Anholt, 1994); and (c) larvae cannibalise both eggs and smaller larvae (e.g., cape anchovy, Brownell, 1985; and many fish species-see review in Pereira, Agostinho, & Winemiller, 2017). Understanding the potential for intra-cohort cannibalism to induce phenological shifts in the optimal breeding schedule can unveil the potential for adaptation under environmental change and explain what otherwise seem to be suboptimal behavior.
We employ this as a basic model by assuming that older individuals are also larger and describe the three different age-dependent cannibalistic modes: egg cannibalised by larvae, larvae cannibalised by elder larvae, and these combined. Our model describes deterministic population dynamics through early life stages and to the age at maturation a m . Transition to the larval stage occurs after the egg development is completed, at age a r . Schematic diagrams of the model are provided in Figure 3.
The dynamics of a population with age a at time t in an age-dependent cannibalistic predation model is thus described by, where μ(a) is the age-specific mortality composed of the natural and cannibalistic mortalities, where μ 1 and μ 2 are the egg and larvae natural mortality rates, respectively, c 1 is the egg cannibalism rate by all larvae, and c 2 is the larval cannibalism rate by elder larvae. Equation (1) is prescribed with the boundary condition of the total number of eggs spawned at time t.
where the total number of egg E = ∫ n(0,t)dt in the breeding season is held fixed (Hastings, 1987;Iwasa et al., 1983). This also provides an opportunity to perform the experiment with the same number of eggs (Hastings, 1987).
The match-mismatch hypothesis normally refers to breeding phenologies adapted to a peak in some key food resource. Here, to keep the analysis simple, we let the abiotic environment drive the hatching success of an egg under the given environment at the day of birth. For example, Reglero et al. (2018) observed that the probability of egg hatching in bluefin tuna varies with temperature and that temperature and therefore egg hatching success in their natural spawning followed a cyclic pattern over the season. The dashed line in Figure 1 includes this and several other temperature-dependent processes, which this represents any abiotic or density-independent seasonal factor that define a breeding window in time.
F I G U R E 2 Schematic diagram of the abiotic environmental seasonal cycle and probability of recruitment success p(t) of two imaginary breeding seasons A and B where the duration of breeding seasons T e are indicated by bidirectional horizontal arrows. Eggs can hatch when environmental condition is above the physiological tolerance limits in an environment (thus, egg hatching success of an egg produced at this time is larger than zero p(t) > 0; the physiological tolerance limits of two hypothetical breeding seasons A and B are represented by E A and E B , respectively, on yaxis). The star icon represents the optimal environmental conditions for an egg to hatch, and without cannibalism, we expect breeding to occur around this time to match this optimal environment with the hatching. With cannibalism, there can be a trade-off between breeding at times of high egg hatching success or earlier, at times with lower risk of cannibalistic predation by older larvae. Therefore, there is a game between parents for optimal breeding time to maximize the survivorship of their own eggs to the adult stage We assume that the number of newly hatched eggs at time t is determined by the eggs spawned at time t − a r , where a r is a fixed egg development time. Moreover, eggs experience the natural mortality rate μ 1 during the egg development period and the egg hatching success is determined by the environmental condition (e.g., temperature or humidity) at time of hatching. Hence, the number of new larval recruits at age a r is, where p(t) is the probability of hatching success at time t. We assume a cyclical environment as in Figure 2, and we describe the hatching probability p(t) over the season as, where nonnegative parameters ∈ 0,1 and ω determine the shape of the cycle, and τ is an arbitrary constant to adjust an environmental peak and the duration of the breeding season. Here, without loss of generality, we set α = 1. The period T p defines the potential breeding period (i.e., when p(t) > 0).
The situation with only egg or larval cannibalism arises by setting c 2 = 0 in Equation (1) and c 1 = 0 in Equation (2), respectively. In our model, when no cannibalism occurs, the optimal time for breeding is independent of the decisions of others and it is a single point in time to match the best environment for egg hatching (Figure 2). determined by the density-independent abiotic conditions. Let us define the egg fitness at breeding time t 0 as the fraction of eggs that survive to the adult stage at time t 0 : When many parents compete with each other to maximize the fitness of their offspring in a noncooperative game-theoretic context or ESS, the equilibrium egg fitness satisfies (Iwasa et al., 1983).
where λ is a positive constant. At this equilibrium, we anticipate that egg phenology in the potential breeding season, T p , shows an emergent egg distribution with a length T e . For convenience, we define the fraction of effective breeding season as T e /T p . The deterministic population dynamics model we use (Equation 1) does not distinguish between single and multiple reproduction events of single individuals in one breeding season, but merely search for an egg distribution that satisfies the condition Equation (6).
As in Iwasa et al. (1983), we assume a constant number of eggs and do not consider the population dynamics of adults. The ESS breeding phenology only represents one single breeding season, since following an ESS in a dynamic population over time involves large computational challenges. To find the optimal breeding schedule under size-dependent cannibalism, we numerically integrate Equation (1) and perform a heuristic algorithm to find the condition satisfying with Equation (6).
See Appendix for more details on the numerical details.

| RE SULTS
The intensity of cannibalism can differ between life stages, and we examined parameter sets representing larvae consuming eggs or smaller larvae or both. We also varied the length of the viable breeding season, and the total number of eggs produced by the parental population. Table 1 is the list for the parameter values examined.
( In general, cannibalism spread the emergent egg distribution over the breeding period (i.e., spawning asymmetry), rather than at a single breeding time as in a noncannibalistic population (Figure 4a-c, top). Also, breeding occurs as soon as the environment allows eggs to hatch (an event that takes place after the egg development time a r ) and most of the breeding takes place earlier than the environmental optimum (stars on x-axis). Hatching earlier is favored as the cannibalism intensity increases, and this trend shifts the distribution earlier in the season, leading to a narrower breeding phenology due to an accumulation of egg density toward the earliest possible breeding. This creates a monotonically decaying breeding distribution. The intensity of cannibalism is also characterized by the egg fitness to the age at maturation, a m (Figure 4a-c, bottom). Cannibalism on eggs causes lower egg fitness, and earlier hatching is favored, and more than for larval cannibalism. The emerging breeding period also varies with the strength of cannibalism, and it is more contracted with higher cannibalism.
We observed similar trends when the environment suitable for breeding is increased to 50% of a year (Figure 4d-f). Counterintuitively, the longer potential hatching periods tend to result in shorter breeding periods. An exception is when cannibalism is high and only on larvae ( Figure 4f). In this case, the distribution of the breeding period is wider when the environmentally suitable breeding window is longer, giving a bimodal breeding distribution with one peak close to the abiotic environmental optimum and one peak early in the season, a signal of a breeding cycle within a single breeding period. Average values and fractions of effective breeding seasons are provided in Table 2, and the qualitative discussions above are verified.
Qualitatively similar patterns appear with different total numbers of eggs, E (see Figure A1 for E = 10 3 and 5 for E = 10 5 , respectively, in Appendix). The exception is again when the cannibalistic intensity is high, and it is only imposed on larvae, and if the total number of eggs is ten times larger (N = 10 5 ; Figure A2c,f), then the smaller egg fitness yields a wider emergent breeding period.

| D ISCUSS I ON
Cannibalism on egg, larvae, or both from larvae hatched earlier in the season leads to intraspecific competition over the optimal breeding period, and apparently, suboptimal breeding phenologies relative to some abiotic optimum may emerge. Firstly, cannibalism causes a breeding asymmetry rather than a simple match to the best seasonal environment. Emergent breeding phenologies are the outcome of a game-theoretic intraspecific competition between parents where each individual makes a decision based on decisions of the others. It differs from other mechanisms inducing hatching asymmetry, such as conservative bet-hedging (Simons & Johnston, 2003) and diversified bet-hedging strategies (Schindler, Armstrong, & Reed, 2015) where the main goal is adaptation to variable environments, not involving the choice of con-specifics (Starrfelt & Kokko, 2012). Secondly, the model predicts a shift of breeding toward earlier and less favorable conditions for the progeny if cannibalism is an important source of mortality. These two properties agree with observations of the spawning phenology of Atlantic bluefin tuna in the Mediterranean Sea (Reglero et al., 2018), where size-dependent larval intra-cohort cannibalism is known to occur (Reglero et al., 2011). Also, De Block et al. (2005) observed early and log-normally distributed breeding phenologies under egg cannibalism in damselfly populations.   TA B L E 2 Mean value and the fraction of effective breeding season, T e /T p , of each curve in Figure 4 (mean, T e /T p ) theory should include dynamics of the full life cycle over multiple generations (e.g., Varpe, Jørgensen, Tarling, & Fiksen, 2009). In addition, a variable environment can lead to bet-hedging reproductive strategies (Schindler et al., 2015;Starrfelt & Kokko, 2012), and under such conditions, the realized phenology may deviate from our deterministic predictions. To assess these factors, we need models of greater complexity, and to connect game-theoretic and bet-hedging concepts.
Our minimal approach, assuming cannibalism causes only an additional mortality, predicts an early shift in breeding phenology. It does not consider the benefit of energy gain by cannibalism, which may further promote an even earlier shift of breeding phenology. A physiologically structured population model is one way to explicitly include the energy gain from cannibalism (Huss et al., 2010), although the model requires greater model complexity and larger number of parameters.
Our model can also help us understand context-dependent trajectories of species assemblage and priority effects (De Meester et al., 2016;Fukami, 2015). Species emerge in seasons in a constant game with other species, and the interactions are size-and time-dependent. Extensions of our modeling approach can include size-dependent competitive interactions between several species and thus become a tool to integrate evolutionary dynamics of phenologies into analysis of environmental change.
Phenological response to climate change differs in a complex manner between species, and it can lead to a trophic mismatch (Both, Van Asch, Bijlsma, Van Den Burg, & Visser, 2009;Edwards & Richardson, 2004). Yang and Rudolf (2010) emphasized the importance of understanding the interaction between phenology and stage-structured species interaction to predict the response to climatic change. Our model is only a starting point, but demonstrates the role that size-dependent intraspecific interactions can have in forming life histories and breeding phenology.

ACK N OWLED G M ENTS
We thank M. H. Cortez and T. Fung for reading a draft and providing thoughtful comments. We also thank two anonymous reviewers for their constructive comments.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

E TH I C A L A PPROVA L
The research complies with all national and international ethical requirements.

DATA AVA I L A B I L I T Y S TAT E M E N T
The manuscript does not contain original data. Here, we present the numerical algorithm that heuristically searches the distribution of egg production that satisfies the condition in Equation (6). The algorithm is based on the following simple steps:

O RCI D
1. Prepare an arbitrary initial egg distribution where the total eggs are E over the potential breeding period b t i Δt = n t i , where b t i Δt is the number of eggs between the time t i and t i + Δt and ∫ b(t)dt = E holds. (1) with the initial egg distribution and calculate fitness at each breeding schedule ϕ(t).

Perform a numerical integration of Equation
3. Find the maximum ϕ max and minimum ϕ min fitnesses among the possible breeding schedules and replace the number of eggs dE from the time with the lowest fitness to the highest fitness.
4. Return to step 2 until the relative difference between ϕ max and ϕ min becomes smaller than the predetermined threshold To avoid a cyclic egg replacement between two and more breeding times, we shuffle the egg number dE between two randomly chosen breeding times for every 10 4 iterations. Moreover, we calculate the relative difference of the average relative fitness difference defined in step 4 for every 300 iterations over 100 iterations, and decrement dE by 5% when this relative difference becomes smaller than the TA B L E A 1 Mean value and the fraction of effective breeding season, T e /T p , of each curve in Figure A1 (mean, T e /T p )

TA B L E A 2
Mean value and the fraction of effective breeding season, T e /T p , of each curve in Figure A2 (mean, T e /T p )   Values used in this algorithm (e.g., predetermined thresholds) are heuristically determined. We used a uniform initial egg distribution.
We observed the order of iterations to obtain convergence spans 10 4 -10 6 depending on the cannibalistic parameter used in the main text. We choose this algorithm because we observe that our objective function i max − i min i max largely fluctuates over iterations (i.e., not monotonically decreases with i), and therefore, ordinary optimization algorithms cannot be directly applied. It is also noted that we truncated the fitness value where an egg load is smaller than 0.5 and define the interval of population n t (a,t) = max n t (a,t) ,0 to facilitate the convergence of our heuristic algorithm.

F I G U R E A 2
The distributions of the normalized egg density (top of each panel: left-vertical axis) under three cannibalism target (i) egg cannibalised by larvae (Egg); (ii) larvae cannibalised by larger larvae (Larvae); and (iii) larvae cannibalise both eggs and smaller larvae (Egg + larvae), and egg recruitment success in relation to the environmental condition (right-vertical axis). Bottom of each panel is egg fitness. These are examined under three cannibalistic scenarios: low (left; c 1 = c 2 = 0.0001); intermediate (center; c 1 = c 2 = 0.001); and high (right; c 1 = c 2 = 0.01). The stars on the x-axis on the top of each panel represent the optimal spawning time under the no cannibalistic predation. The potential breeding period is T p = 0.25 (top two rows) and T p = 0.5 (bottom two rows). The number of eggs is E = 10 5 . Other parameter values are shown in