Thermal pace-of-life strategies improve phenological predictions in ectotherms

Phenological variability among populations is widespread in nature. A few predictive phenological models integrate intrapopulational variability, but none has ever explored the individual strategies potentially occurring within a population. The “pace-of-life” syndrome accounts for such individual strategies, but has yet to be explored under a phenological context. Here we integrated, for the first time, the slow-fast thermal strategies stemming from the “pace-of-life” into a mechanistic predictive framework. We obtained 4619 phenological observations of an important crop pest in the Bolivian Andes by individually following 840 individuals under five rearing temperatures and across nine life stages. The model calibrated with the observed individual “pace-of-life” strategies showed a higher accuracy in phenological predictions than when accounting for intrapopulational variability alone. We further explored our framework with generated data and suggest that ectotherm species with a high number of life stages and with slow and/or fast individuals should exhibit a greater variance of populational phenology, resulting in a potentially longer time window of interaction with other species. We believe that the “pace-of-life” framework is a promising approach to improve phenological prediction across a wide array of species.

Over the last decades, climate change has induced an advance in phenological events for many plant and animal species [1][2][3] . Accurately predicting phenological changes is challenging because of the interactions between numerous causal factors driving phenological events [4][5][6] . Among these factors, phenological responses between individuals can account for a significant amount of variation in phenology within populations, even though this variation has seldom been studied 7 . A non-exhaustive review of existing literature (Supplementary Table S1) reveals that phenological differences between individuals are a widespread phenomenon across taxa, and are in part explained by the mosaic nature of the environment, with each individual experiencing a unique set of environmental conditions 8,9 . Inter-individual differences in phenology also occur under controlled conditions, suggesting that physiological drivers are also involved. In particular, inter-individual differences in development affect phenology because the duration of a life stage at the population level depends on (i) the inter-individual differences in synchrony of the life stage, and (ii) the inter-individual differences in time needed to complete the stage 4,10 . Even though some phenological models integrate intrapopulational variation in development 11,12 , none has ever explored the individual strategies potentially occurring within a population.
The "pace-of-life" (POL) syndrome postulates that a population can be divided into subpopulations with "fast" and "slow" individuals reflecting their life-history strategies and adaptations to different environmental conditions 13 . Each pace displays several correlated traits, such as energy expenditure, growth rate or lifespan 14 . Such individual strategies in performance rate occurs in natural population and has been reported for several species of trees 15 , birds 16,17 , mammals [18][19][20] and invertebrates 21,22 . However, POL strategies have seldom been applied to development studies (but see 22 ) as this requires the individual following of several populations over their entire life cycle. Life cycles of most organisms include different periods (referred as life stages) throughout their ontogeny, with stage-specific form, lifestyle, reproductive capacity, or physiology 10,23,24 . Thus, differences in development rate not only occur between individuals, but also between life stages within an individual's life cycle. The POL SCIENtIFIC REPORTS | (2018) 8:15891 | DOI: 10.1038/s41598-018-34274-1 syndrome could be divided across the several life stages, with faster and slower individuals for a specific life stage. However, it is unclear whether the development strategies are consistent across life stages 22 . Surprisingly, little attention has been paid to the incorporation of POL strategies in phenological predictions, despite a novel focus on POL strategies in other performances [15][16][17][18][19][20][21] .
Inter-individual differences in thermal performances have been suggested as a source of the POL syndrome, through their effects on metabolism 25 . Ectothermic organisms are ideal models to study this relationship as environmental temperature constrains several physiological performances including development that acts as a major determinant of the individual's phenology 26,27 . Moreover, each life stage exhibits different ranges of performance responses to temperature (usually characterized by thermal performance curves) that reflects its adaptation to a specific habitat 28 . Even when experiencing constant temperatures, ectotherms show variability in emergence dates, owing to differences in development rate across individuals and life stages. However, when submitted to fluctuating temperatures, complex interactions can arise between the physiological differences and the temperatures experienced, with potential lags and accelerations in development time between individuals 10 . Lag between individuals implies that different temperatures can be experienced during a specific development event, which has long-lasting consequences on future development and metabolism 29,30 .
We hypothesize that taking into account thermal POL strategies in development rate models would improve phenological predictions in ectotherms. To address this issue our objectives are to: (i) determine the occurrence of thermal POL strategies in development, (ii) develop and assess a new framework accounting for thermal POL strategies in phenological models, and (iii) explore the effect of population slow-fast strategies composition on phenological predictions.
We developed three individual-based models following a growing complexity, and compared them using real development data from Copitarsia incommoda. Model 1 employs the commonly-used thermal performance curve that estimates a continuous mean performance rate for the whole population 31 , resulting in identical phenology across individuals (Fig. 1A2). Model 2A overcomes this issue by capturing a range of development rates naturally occurring within a population (Fig. 1B). Finally, Model 2B accounts for individual slow-fast thermal strategies ( Fig. 1C; see Methods for details about the models).

Results
Models' predictions comparison: predicting the phenology of the quinoa moth. The proportion of individuals with slow and fast strategies within each population of C. incommoda was similar, and accounted for 37 to 59% of the population (18 to 29% for slow strategies and from 17 to 30% for fast strategies; Fig. 2). By accounting for these strategies, the predictive ability of Model 2B outperformed Model 2A in 89% of cases (  Table S3).
Phenological variance and windows of interaction. We found that the phenological variance (i.e. the distribution of emergence dates) for the last life stage increased with the proportion of slow and fast strategies occurring within the populations, and with the number of total life stages (Fig. 4). A population from a species with four life stages composed of 10% of fast and 10% of slow individual strategies showed a normalized phenological variance of 0.1, whereas the normalized variance widened up to threefold when half of the individuals follows a fast strategy and the other half a slow strategy. A species with eight life stages experienced a similar threefold increase between 10-10% and 50-50% slow-fast strategies but with a wider overall variance (Fig. 4, Supplementary Fig. S3).

Discussion
By independently following the individuals of C. incommoda during their whole life cycle, our study revealed the occurrence of pace-of-life (POL) strategies in the development of an ectothermic species. The slow and fast individuals accounted for a large part of the population (between 59 and 36%), with a roughly equal repartition between slow and fast (Fig. 2). As reported in other studies 32,33 , our population was not composed solely of strict slow and fast individuals but also included individuals showing an intermediate strategy with a phenotype close to the average of the population. This was due to individuals shifting between slow and fast strategies between life stages (Fig. 2). The lack of repeatability for performances (e.g. metabolism, behavior) across life stages has been observed in species of mammals 34 , birds 35 , reptiles 36 , amphibians 37 and arthropods 38 . The individual POL strategies are more likely to switch between life stages with strong physiological differences, such as life stages before and after metamorphosis in amphibians or arthropods 37,38 . However, in the case of C. incommoda, we could not find any consistent patterns across individuals (Fig. 2). This result could be linked to the fact that the individual life-history has a strong impact on development compared to other performances that happen at a shorter timespan (e.g. locomotion). Moreover, because the life stage progress is linear, the inter-individual differences occurring at one life stage affect the outcome of the next life stage 10 . As the rearing temperatures were slightly fluctuating, differences in the temperatures experienced during early-life may have further affected the development during subsequent life stages. These peculiarities of development compared to other performances further support the development of thermal performance probability models in phenology studies.
A high predictive ability under several thermal contexts was shown in Models 2A and 2B (Fig. 3). Nevertheless, Model 2B performed better by integrating the slow-fast strategies stemming from the POL than Model 2A that accounted for intra-populational variability in development rates alone (Fig. 3 models account for inter-individual differences in development to predict the range of phenological responses observed in nature 12,39,40 . To the best of our knowledge, this study is the first published to integrate the POL into an individual-based phenological modelling framework. The POL syndrome has been widely used to link an organism's behavior with its metabolism 20,41,42 . Recently, is has been argued that this syndrome arises from lower level processes such as thermal physiology 25 and that it is relevant at the individual level 14 . We showed that the inter-individual differences and the POL strategies are also relevant to predict phenology. The accurate quantitative predictions generated by our model strongly supports the importance of POL strategies as a driver of phenology 43 . Our simulations suggest that the proportion of slow and fast individuals within a population of virtual ectotherms affects the phenology (c.f. distribution of individual emergence dates) of the population (Fig. 4). Larger proportions of slow or fast individuals increase the phenological variance, and the combination of slow and fast individuals further widen the variance (Fig. 4). The overall pattern arising from the simulation highlights a complementarity effect between the slow and fast strategies, which may have important ecological consequences. Indeed, the variance of a populational phenology delineates the potential interaction window with another species, which often change in strength and type across ontogeny 4 . Therefore, a phenological mismatch between ontogenies can result in weaker or shorter interactions 44 , or even change the type of interactions between species (e.g. predation to competition) 45 . In this context, we suggest that populations with a higher proportion of slow and fast individuals are expected to be more resilient to temporal mismatches between species, as their phenological window is wider. Our results also suggest that, for a given proportion of slow-fast individual, the number of life stages increases phenological variance (Fig. 4). As the number of life stages in ectothermic organisms shows a high variability across taxa 23 , temporal mismatch between species is likely to be stronger for species with fewer  49 . Spatial thermal heterogeneity due to biotic and abiotic factors also acts as a driver of thermal fluctuations experienced by small ectotherms moving across a mosaic of microclimates 50 . Therefore ectothermic animals can buffer thermal fluctuations by remaining in microhabitats near their optimal temperature 51,52 , or through behavioural thermoregulation 53 . A step further into the understanding and prediction of the inter-individual variations in phenology in natural habitat should integrate both these behavioural and spatial fluctuations in temperature.

Methods
Phenological models. Phenological predictions for a population of ectotherms (such as reptiles, amphibians or arthropods; Fig. 1) can be modelled by (i) determining the development response of the population to temperature, and (ii) accumulating the development rate over a temperature time-series 54 . A standard approach to studying the relationship between development and temperature is carried out by measuring the development rate across a range of environmental temperatures under controlled conditions (i.e. rearing experiments). The performance rate typically rises slowly with increasing temperatures until it reaches a maximum performance and then drops more or less abruptly (Fig. 1A). The accumulation process for each individual consists of (i) following a temperature time series with a specific time step (Fig. 1A1,B1), (ii) retrieving the development rate  corresponding to the temperature at one discrete-time frame (triangles, squares and circles in Fig. 1A,A1), (iii) multiplying the development rate by the time step, (iv) accumulating the development rate until the next life stage is reached (e.g. juvenile lizard, tadpole and insect larva in Fig. 1A2), and (v) recording the time elapsed during the whole accumulation process (Fig. 1A2,B2,B3). The modelling steps are then repeated along the life stages (colours in Fig. 1) and individuals (to obtain the distributions of phenology in Fig. 1A2,B2,B3). Starting with this general framework, we propose increasingly complex models by deepening the underlying hypotheses behind the thermal characterization and the development accumulation processes. Our models account for phenological events during growth, but set aside the seasonal timing of such events.
Model 1 -The commonly used thermal performance curve. Model 1 employs the commonly-used thermal performance curve (TPC) that allows estimating a continuous mean performance rate for the whole population 31 . The fitting process can be achieved by using simple or more complex models 31 . Once the TPC for development has been characterized, a development rate can be determined at any temperature within the thermal limits of the species (Fig. 1A1,A2). However, because the TPC offers a mean development rate response (Fig. 1A1), the dates of emergence for every life stages are identical across individuals (Fig. 1A2).
Model 2A -Thermal performance probability. Model 2A is motivated by the fact that development occurs at longer time scales than other performances (e.g. locomotion). Indeed, for a given life stage, development rate cannot be measured at different temperatures on the same individual, which impedes the application of models quantifying within versus between individuals variances 55,56 . We chose a modelling approach easily integrated into the well-established performance accumulation models from TPC 27 . Instead of fitting a curve close to the mean of population measurements at various discrete temperatures, we propose to fit a mixture distribution of development responses at each discrete temperature 11 ( Supplementary Fig. S1). Finite mixture distributions are commonly used to identify sub-populations within a population 57 . These models assume that the sample observations randomly arise from two or more distributions with certain probabilities. Suppose that = ... X X X ( , , ) n 1 is a random sample of size n from the mixed population, with a density function of X i given by f(X i ). In a generic case, X is assumed to have arisen from a mixture distribution with two components (K = 2) following normal distributions. Therefore, the density of X i is given by: Where λ k is the weight, μ k is the mean, and σ k 2 the variance of normal component k. We used an expectationmaximization (EM) algorithm to identify the mixture distribution parameters at the discrete known temperatures. The EM algorithm is decomposed in two steps, the expectation (E) step based on Baye's theorem and the maximization (M) step (normalmixEM function from the mixtools R package) 57 . Several studies suggested that the distribution of development rates within a population follows a Weibull 40 or a log-normal distribution 11 . However, following our pace-of-life hypothesis, we decided to fit a mixture of two normal distributions (i.e. slow and fast individuals). We found higher log likelihood scores for mixture distributions (i.e. bimodal distributions) than for Weibull and lognormal distributions in 42 out of 45 cases (Supplementary Table S2), therefore supporting our hypothesis.
To determine development rate distributions between the measured temperatures, we interpolated the normal distributions' parameters (mean, standard deviation, and weight) to obtain a continuous response over the whole range of temperatures ( Supplementary Fig. S1). Each parameter (p) from the mixture distributions was interpolated between two known development distributions at consecutive temperatures with the following expression: Where T 1 is the lower temperature at which parameter p T 1 is known, T 2 the higher temperature limit for the interpolation at which parameter p T 1 is known, and T x corresponds to the temperature at which the parameter (p(T x )) has to be determined. The parameters interpolation is then repeated for every parameter (λ k , μ k , and σ k 2 ) and for the two normal components of the mixture distribution. We assumed that a linear interpolation would offer a sufficient fit if the number of rearing temperatures is high and evenly distributed across the thermal range. However, if the number of temperatures is low or does not incorporate the thermal limits, we advise using a non-linear interpolation based on one of the TPC models available, allowing to determine the thermal limits.
We named the resulting probability surface the Thermal Performance Probability (TPP), which is depicted graphically as a heat map with its area delimiting the response space, and the colour density representing the probability of the response (Fig. 1B, Supplementary Fig. S1). During the development accumulation process, more than one development response is possible at a given temperature, and the determination of this development rate is performed by drawing one value from the mixture distribution. This approach allows capturing a range of theoretical responses naturally occurring within a population (Fig. 1B). It results in a distribution of phenology for each life stage overlapping each other, with potential accelerations or buffering effects (Fig. 1B1,B2,B3).

Model 2B -Thermal performance probability with pace-of-life strategies. A limitation of Model
2A is that one individual can exhibit development rates that vary from one extreme to another at two consecutive time frames within the same life stage, which may be unrealistic in natural populations, as suggested by the POL syndrome. To address this limitation, we implemented individual slow-fast strategies in Model 2B by dividing

Study case.
To compare and illustrate the three models, we used development data from the quinoa moth (Copitarsia incommoda), and predicted its phenology distributions for each life stage. C. incommoda is an important pest of the quinoa crop (Chenopodium quinoa) and shows several characteristics representative of many species across the world. First, it thrives in a thermal environment with strong daily fluctuations (pers. obs.). Because a trade-off exists between adaptations at high and low temperatures, a high variability in performance responses to temperature is expected within and among the populations 27 . Second, the quinoa moth has numerous life stages (egg, six larval instars, pupa and adult) and each life stage has a particular thermal environment, which may trigger variable development rates across life stages 58 . Third, the quinoa moth lives in an environment with virtually unlimited food supply, which lowers the importance of this confounding factor on development. We reared each C. incommoda individual from each population in separate containers with identical thermal conditions to daily follow its development. By doing so, we were able to identify potential slow-fast strategies occurring among individuals and calibrated the models accordingly. We then compared models' predictions with the observed phenology using a cross-validation approach.
To obtain data about the relationship between temperature and development, 960 eggs of C. incommoda were distributed among 8 rearing units with different temperatures spanning the range of temperatures in natural conditions (mean temperatures of 5.1, 12.6, 18.1, 20.2, 24.8, 30.0, 33.2, and 34.6 °C and standard deviations of 1.5, 1.3, 1.8, 2.2, 1.6, 1.0, 0.9, respectively, with 60 ± 5% relative humidity). We chose to create temperatures slightly fluctuating around a mean with low amplitude instead of strictly constant temperatures in order to stimulate the differences in development accumulation between individuals, without facing the issues of highly fluctuating temperatures 48 . Each rearing unit was checked once a day following a regular schedule, in which process life stage of each individual was individually followed (total of 4619 phenological measurements). Temperature in each rearing unit was measured every 30 minutes using temperature loggers (HOBO TidbiT v2 UTBI-001 from Onset) in order to feed the model with actual temperature time-series. Detailed information about the rearing experiments can be found in Rebaudo et al. 58 . The experiments complied with the Association for the Study of Animal Behaviour guidelines for the use of animals in research 59 .
To construct the models and evaluate their predictions, the development rate dataset for each temperature was split into two datasets following the cross-validation approach 7 : a calibration dataset (70% of the observations) and an evaluation dataset (30%). The three models were fitted on the development rates of the eight populations from the rearing units at eight mean temperatures. However, we used only the five temperatures where egg mortality was lower than 50% in the subsequent modelling steps (12.6, 18.1, 20.2, 24.8, 30.0 °C). We fitted Model 1 for each life stage using the Lactin-1 model 60 to obtain the predicted phenologies for each life stage and to compare it to the observed one. All TPCs were computed using the devRate R package 61,62 . To construct Model 2B, we determined the empirical contribution of fast, slow, and intermediate individuals within each population by splitting the development range in two considering the median development value. We then followed each individual from egg until the sixth (last) larval instar, and recorded in which category (i.e. slow, or fast or intermediate relative to the median) it fell at each life stage. We decided to attribute a slow or fast strategy to an individual considering its strategy consistency across life stages. To do so, we determined the proportion of fast and slow individuals in each population by separating i) the fast individuals that stayed in the upper category for at least 5 out of 7 life stages (population reared at 13 °C = 30%, 18 °C = 17%, 21 °C = 19%, 26 °C = 20% and 30 °C = 22%) and ii) the slow individuals that stayed in the lower category for at least 5 out of 7 life stages (population reared at 13 °C = 29%, 18 °C = 19%, 21 °C = 18%, 26 °C = 24% and 30 °C = 22%). We focused on the larval stages because of the high mortality during the subsequent life stages. The three models were run along five independent temperature time-series (time-step of 24 h). The goodness-of-fit between predicted and observed phenology distributions were evaluated differently for Model 1 and Models 2 according to their outputs. Model 1 was evaluated by computing the Root Mean Squared Error between predicted and observed data as it does not account for variance in phenological outputs. Model 2A and 2B were evaluated by computing the overlap that accounts for the variance in the output. Measuring the overlap goes beyond comparing means and variances. An overlap of 1 means complete overlap between two distributions (i.e. perfect prediction), while an overlap value of 0 indicates an absence of overlap. The overlap estimator Dhat1 was used because most of our sample sizes were lower than 50 63 . The overlap were estimated for Model 2A and 2B using the overlapEst() function from the overlap R package 64 . All analyses were compiled with R 3.3.3.
Exploring the thermal performance probability framework. To explore how Model 2B parameters affect phenology patterns, we generated development data for a virtual species with a varying number of life stages (4 and 8) and slow and fast proportion in the population (0, 10, 20, 30, 40, and 50% for each strategy). The generated population of 100 individuals had an inter-individual variability in development rate with an identical TPP across life stages. The model was performed on generated constant temperature time-series. We then assessed the combined consequences of these parameter values on the variance of the phenology distribution for the last life stage of the whole population, based on 30 runs for each combination of slow and fast strategies.

Data Availability
All data and computer source code used to calibrate and run the models are accessible with the https://doi. org/10.5281/zenodo.1421887.