Optimal reproductive phenology under size-dependent cannibalism

Cannibalism is a size-mediated priority effect, and introduces a tradeoff between breeding at the optimal time in the season and increased predation risk from offspring born earlier. This game-theoretic situation among parents may drive adaptive reproductive phenology. However, it is not straightforward to quantify how cannibalism affects egg fitness or to distinguish emergent breeding phenology from alternative adaptive drivers. Here, we employ an age-structured game-theoretic mathematical model to find evolutionary stable breeding phenologies. We predict how size-dependent cannibalism on eggs, larvae, or both change emergent breeding phenology, and find breeding that 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. This suggests that the nature of cannibalism among early life stages is a potential driver of the diversity of reproductive phenologies seen across taxa.


Introduction
The reproductive schedule of species has a tight connection to fitness and ontogenetic development [1], and understanding the interplay between breeding phenology and size-dependent interactions is needed to understand how organisms adapt to a changing climate [2]. Improving our ability to predict phenologies and seasonal structure of species interactions not only facilitates our view of evolutionary traits but also provides significant merit for ecosystem management, conservation, and estimations of the impact of climatic change on species [2][3][4].
The colonization of habitats and arrival times of species in a seasonal environment is important in predator-prey interactions and in forming the structure of communities [5][6][7][8].
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 [5], as organisms emerge small and grow larger while their role as predators and prey shift [9]. An early start can be a large advantage both in competitive and in predator-prey interactions, but are traded against the match with key resources during ontogeny [10].
Cannibalism is a size-mediated priority effect that occurs within species. It is ubiquitous in many groups of animals, including insects, fish, amphibians, birds, and mammals [11,12], and can amount a major part of early-life mortality [13,14]. Intraspecific oophagy is rather common in many egg-lying animals, and newborn individuals are similarly vulnerable to cannibalism due to their limited ability to avoid predation and nutrient rich composition [12]. It also drives population dynamics [15][16][17][18] and reproductive behavior [11,12,19]. If there is a seasonal peak in the environmental growth or survival conditions for eggs and larvae, cannibalism during early life stages introduces a tradeoff 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 tradeoff may induce various adaptive breeding phenologies observed across taxa, for example the spatiotemporal split of larval load in the fire salamander population [20], the surprisingly early and low-temperature spawning period of the Atlantic bluefin tuna [21], or the early hatching of damselfly [19]. 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 non-cooperative game-theoretic approach to find the evolutionary stable strategy (ESS) [22,23] under size-dependent cannibalism. To our knowledge, however, there is no theoretical studies of this process to promote our understanding of size-mediated priority effects.
Here we employ a theoretical approach to investigate how size-dependent cannibalism can shift breeding phenology in a seasonal environmental cycle. The common starting point to understand the timing of egg production 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 [10,24]. 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. Ideally, the critical period concept should include the full future life cycle of the egg at the time it is laid (e.g. [25]), but this require a full life-cycle model.
If cannibalism is sufficiently strong, it can lead to deviations from the expectation of all offspring appearing at the optimal environmental conditions [21]. We define an optimal time for breeding based on the environment (Fig. 1), but add cannibalism among early life stages and game-theoretic decision making of breeding time among parents. The intensity of cannibalism is quantified by survivorship through the egg and larval stages (egg fitness).
Here we intend to achieve general insight into the effect of cannibalism on breeding phenology and, for this purpose, we discuss three general size-dependent cannibalistic modes: (i) egg cannibalised by larvae (e.g., anuran [26] and Tribolium (flour beetle) [17]); (ii) larvae  Figure 1: Schematic diagram of the environmental seasonal cycle and probability of recruitment success of two imaginary species A and B. Eggs can hatch and survive to adult stage when environmental condition is above the physiological tolerance limits (thus, egg hatching success of an egg produced at this time is larger than zero; the physiological tolerance limits of species A and B are represented on y-axis). The sun icon represents the best environmental conditions for an egg to hatch, and we expect breeding to occur around this time when there is no other constraint on the spawning decision making of parents. There can be a tradeoff between high egg hatching success and low risk of cannibalistic predation by older larvae. Therefore, there is a game between parents to find an optimal breeding time to maximize the survivorship of an egg to the adult stage cannibalised by larger larvae (e.g., damselfly [14]); and (iii) larvae cannibalise both eggs and smaller larvae (e.g., cape anchovy [27] and many fish species). Understanding the potential for cannibalism to induce phenological shifts in the optimal breeding schedule can unveil the potential for adaptation to environmental change and explain what otherwise seemingly be non-optimal behavior.

Material and Methods
Population dynamics that incorporate age-dependent cannibalism is often described by the McKendrick-von Foerster equation [15][16][17]. 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 population dynamics through early life stages and to the age at maturation a m . Recruitment occurs at the age at recruitment, a r , after the egg development time t r . Schematic diagrams of the model are provided in Fig. 2.
The dynamics of a population with age a at time t in an age-dependent cannibalistic  Figure 2: Schematic diagram of size-dependent cannibalism in the model. A constant rate of mortality from other sources is also included. We do not explicitly consider the dynamics of adult individuals and, hence, the corresponding age axis is described by the dotted line.
predation model is thus described by where, µ(a) is the age-specific mortality where, µ is the natural mortality rate, c 1 is the egg cannibalism rate by all larvae, and c 2 is the larvae cannibalism rate by elder larvae. Eq. (1) is prescribed with the boundary condition for the 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 [17,22]. This also provide an opportunity to perform the experiment within the same setting [17]. 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 key resource be driving the hatching success of the eggs. For example, Reglero et al. [21] 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 area followed a cyclic pattern over the season. Here, we assume that the number of newly hatched larvae at time t is determined by the eggs spawned at time t − t r , where t r = a r is a fixed egg development time. Moreover, eggs experience the natural mortality rate µ during the egg development period and the egg hatching success is determined by the environmental condition (e.g., temperature or food resource) 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 Fig. 1, 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 merely to adjust an environmental peak and the duration of breeding season. 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 Eq. (1) and c 1 = 0 in Eq. (2), respectively. When no cannibalism occurs, the optimal time for breeding is a single point in time 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 non-cooperative game-theoretic context or ESS, the equilibrium egg fitness satisfies [22] 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 . We define the fraction of effective breeding season as T e /T p . For parsimony, we do not consider the population dynamics of adult individuals. To find the optimal breeding schedule under size-dependent cannibalism, we numerically integrate Eq. (1) and perform a heuristic algorithm to find the condition satisfying with Eq. (5). See Appendix for more numerical details.

Results
The intensity of cannibalism can differ between life stages, and we examined parameter sets mimicking 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. We chose these parameters to facilitate our general understanding of cannibalism's effect on the breeding phenology rather than focusing on a species specific situation. Table 1 is the list for the parameter values examined.
In general, cannibalism spread the emergent egg distribution over the breeding period, rather than at a single breeding time as in a non-cannibalistic population (Fig. 3 (a)-(c), top). Also, breeding occurs as soon as the environment allows eggs to hatch 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 towards 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 ( Fig. 3 (a)-(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 50% of a year (Fig. 3 (d)-(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 (Fig. 3 (c) and (f)). 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 environmental optimum and one peak early in the season. Average values and fractions of effective breeding season are provided in Table 2 and the qualitative discussions above are verified.
Qualitatively similar patterns appear with different total numbers of eggs, E (see Figs. A.1 for E = 10 3 and A.2 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 ; Figs. A.2 (c) and (f)). Then, the smaller egg fitness yields a wider emergent breeding period.  Table 1.

Discussion
Cannibalism on egg, larvae, or both from larvae hatched earlier in the season lead to intraspecific competition over the optimal breeding period, and apparently suboptimal breeding phenologies may emerge. Firstly, cannibalism causes a distribution of the breeding period rather than a simple match to the maximum seasonal resource level. Emergent breeding phenologies are the outcome of a game-theoretic intraspecific competition between parents, and it differs from other mechanisms inducing hatching asymmetry reported previously, such as conservative bet-hedging [28] and diversified bet-hedging strategies [29]. Secondly, the model predicts a shift of hatching towards earlier and less favourable conditions for the progeny if cannibalism is an important source of mortality. These two properties agree with the observation of the spawning phenology of Atlantic bluefin tuna in the Mediterranean Sea [21], where size-dependent larval intra-cohort cannibalism is known to occur [30]. Also, DeBlock et al. [19] observed early and log-normally distributed breeding phenologies under egg cannibalism in damselfly populations. Emergent breeding phenology is sensitive to which life-stage is most susceptible to cannibalism. Our model predicted that egg cannibalism tends to induce an earlier shift than larval cannibalism, since later hatching causes greater cannibalism by the larvae without self-regulation mechanisms due to the lack of intraspecific competition. Hence, if eggs are vulnerable to intra-cohort cannibalism then total predation from cannibals are higher and the drive towards earlier spawning is stronger. For cannibalism acting on larvae, on the other hand, the cannibals are themselves thinned out by cannibalism. This limits the magnitude of cannibalistic mortality among smaller larvae, and allows parents to breed later and under a better environmental conditions. A longer potential breeding period (T p = 0.5) gave a lower ratio of emergent:potential breeding season (T e /T p ) than the shorter potential breeding period (T p = 0.25). In a shorter potential breeding period, the tradeoff between the better seasonal environment and reduced predation risk is weaker, and provides a higher benefit to individuals breeding later in the season than the case of a longer potential breeding period, leading to less contracted breeding phenologies. The magnitude of change in the ratio again depends on cannibalism rate.
Also, depending on the intensity of cannibalism over the ontogeny cannibalism can induce surprisingly diverse patterns of emergent breeding phenologies, including a symmetric bell-shaped distribution with and without cut off at the earliest possible breeding time or monotonically decreasing egg production curves (e.g., Fig. 3 (a) and (c), top). The latter distribution appears in response to intensified cannibalism, and as a rule of thumb, intensifying cannibalism causes a shift of the egg distribution toward an earlier time and it leads to a cut-off of the distribution at the earliest possible breeding time (e.g., Fig. 3 (d)-(f),  top). Even a bimodal pattern (e.g., Fig. 3 (f), top) or more complex phenologies (Figs. A.2 (c) and (f), top) can occur if cannibalism is only on larvae by elder individuals and its effect on survival is high. These diverse patterns were previously reported in multiple taxa. The spawning phenology of Atlantic bluefin tuna resembles a bell-shaped curve [21], while a monotonically decreasing pattern was observed in the burying beetle population under the filial cannibalism [31]. In the context of an optimal breeding schedule of male butterflies to facilitate mating success, Iwasa et al. [22] predicted a bell-shaped pattern with a cut-off as the optimal schedule.
We have used a simple eco-evolutionary model to analyse the effect of cannibalism on breeding phenology. A more comprehensive understanding of cannibalism can be made by explicitly modelling breeding behavior in matured individuals, perhaps for more specific case studies. Provided with a good understanding of the lifecycle and the cannibalistic interactions of concerned species, a wide variety of tactical models can be constructed based on our general model, by, e.g., introducing allometric relationships in fecundity, growth and maturation, or a more mechanistic formulation of the process of cannibalism.
The use of eco-evolutionary models can also help us understand context-dependent trajectories of species assemblage and priority effects [6,32]. Species emerge in seasons in a constant game with other species, and the interactions are size-and time-dependent. Extensions of our modelling 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 [33,34]. Yang and Rudolf [2] emphasized the importance of understanding the interaction between penology 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.
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 predetermine threshold To avoid a cyclic egg replacement between two or more breeding times, we shuffle the egg number dE between two randomly chosen breeding times for every 10000 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 predetermined threshold is the average relative fitness difference where its sampling starts at iteration time i.
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 largely fluctuates over iterations (i.e., not monotonically decreases with i) and, therefore, ordinary optimization algorithms can not 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. 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 3 . Other parameter values are shown in Table 1. 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 Table 1.