Modelling diapause termination and phenology of the Japanese beetle, Popillia japonica

We deve﻿loped a mechanistic, stage-structured model simulating the phenology of Popillia japonica. The model simulates the influence of soil temperature on the larval diapause termination and on the development rate function of post-overwintering larvae and pupae. Model parameters are estimated based on literature evidence for pupae development and on a parameterisation process that allows estimating parameters for larval diapause termination and for the development rate function (and the related uncertainty) of post-overwintering larvae. Data used for model parameterisation and validation refer to time-series adult trap catches collected during the P. japonica monitoring programme performed by the Phytosanitary Service of Lombardy Region within the infested area in Lombardy (Italy) from 2015 to 2019. A total of 12 randomly selected locations are used to estimate biologically realistic model parameters (parameterisation dataset). We applied a Jackknife nonparametric resampling procedure on the parameterisation dataset to quantify uncertainty associated with parameters’ estimates. Parameterised model is then validated on time-series adult trap catches data referring to a different set of 12 randomly selected locations (validation dataset) surveyed in Lombardy. The model successfully predicted the beginning of adult emergence and the overall curve of adult emergence in the validation dataset. The model presented can support the definition of the best timing for the implementation of monitoring and control activities for the local and the area-wide management of P. japonica.


Introduction
The Japanese beetle (Popillia japonica) is an invasive and highly polyphagous pest belonging to the family of Scarabaeidae. The species has been reported on more than 300 Communicated by Jon Sweeney.
The species originated from North-Eastern Asia where it is native in Northern China and in Japan. In 1916, the species was introduced in North America (Fleming 1976) where it became invasive and a serious pest for turf grass (Potter and Held 2002). In the 1970s, the species was found in Europe in Terceira island (Portugal) where it was able to establish and spread on other Azorean islands (Vieira 2008). In 2014, the species was first found in mainland Europe, in Italy in the Ticino Valley along the border between Lombardy and Piedmont Regions (Pavesi 2014). The species was also reported in Switzerland since 2017 (EPPO 2017(EPPO , 2020. Various aspects of the species life-history are influenced by temperature, including the emergence of adult individuals, mating, oviposition and larval development (Fleming 1972). The species is univoltine in the vast majority of its area of distribution, except in the northernmost areas where individuals might take up to 2 years to complete the life cycle (Clausen et al. 1927;Fleming 1972;Vittum 1986). Depending on temperature, adults emerge in summer from May to July (Fleming 1972) and feed on the flowers and fruits of many host plants (Metcalf and Flint 1962;Potter and Held 2002;Klein 2008). Mated females lay from 40 to 60 eggs (Campbell et al. 1989) 10 cm deep in soil. Eggs hatch in 10-14 days, and newly emerged larvae start feeding on the root system of crops and weeds usually in the upper 7.5 cm of soil (Metcalf and Flint 1962). There are three larval instars. The first instar develops in 2-3 weeks and the second instar in 3-4 weeks. Approaching late autumn, second and third instar larvae start moving deeper into the soil (at depths greater than 15 cm) to overwinter and avoid freezing temperature (Potter and Held 2002). In the cold season, larvae enter diapause with a complete inhibition of their development (Ludwig 1928(Ludwig , 1932(Ludwig , 1953. During spring, as soil temperature increases, larvae start migrating to shallower depths, complete their development and pupate. Depending on temperature, pupation might take 1-3 weeks (Ludwig 1932(Ludwig , 1939Fleming 1972;Potter and Held 2002. Currently, the management of P. japonica in Italy is under the official control of the National Plant Protection Organization and involves the Regional Phytosanitary Services, which seek to contain pest populations and prevent the spread of the species. Although control measures were immediate, the species is considered not technically eradicable due to its biological and behavioural characteristics. The management of P. japonica populations might benefit from the development and use of modelling tools supporting rational decision-making for both well established (Gilioli et al. 2016;Rossi et al. 2019) and spreading populations (Gilioli et al. 2013;Ponti et al. 2015;Zhu et al. 2017).
Phenological models are among the most useful tools for predicting the emergence of susceptible life stages (e.g. pre-and post-diapause larvae and adults), thus allowing to schedule monitoring and control actions accordingly (Samietz et al. 2007;Pasquali et al. 2019). Mechanistic models are able to simulate the life-history strategies of the species and how they are influenced by relevant abiotic (e.g. temperature, relative humidity) and biotic (e.g. interspecific competition, plant phenology) drivers (Régnière et al. 1981;Gilioli et al. 2021). Therefore, they represent valuable tools for predicting the phenology of P. japonica.
In this work, we present a temperature-driven, mechanistic model predicting larval diapause termination, postdiapausing larval and pupal development, and emergence of the adult stage of P. japonica. The model accounts for the realistic representation of the role of soil temperature influencing the life-history strategies of the species (Ludwig 1928(Ludwig , 1930Fleming 1972;Gutierrez 1996;Régnière et al. 2012). Data used for model parameterisation and validation were collected by the Regional Phytosanitary Service of the Lombardy Region (Italy) within the infested area in Lombardy, from 2015 to 2019. Model performance in simulating the diapause termination period and in predicting the phenology of P. japonica has been evaluated in terms of model's capacity to describe both the overall shape of the observed adult emergence and the beginning of the first flights calculated at the 2nd, 5th and 10th percentile of adult emergence distribution. These two different aspects related to the phenology of the species are of fundamental importance for supporting the management of the species. Knowing in advance the period of first emergence can be useful for the timely implementation of monitoring actions. Knowing in advance the overall adult emergence pattern allows for the description of the entire adult emergence period and ultimately to estimate the periods of major increase in population abundance associated with the species' impacts. This knowledge can be used for the timely implementation of control actions, increasing the efficacy of the treatment and preventing population growth.

Conceptual model and model assumptions
The model simulates the phenology of post-diapausing larvae and adults of P. japonica, the most relevant part of the life cycle in relation to the management of this pest. The model includes two sub-models: (1) a sub-model (M1) simulating the termination of the diapause period of the larval stage as function of soil temperature triggering the development of post-overwintering larvae and (2) a submodel (M2) simulating the development of post-overwintering larvae and pupae, and the emergence of the adult stage as function of soil temperature. We assume that larvae enter diapause at the beginning of the third instar (Potter and Held 2002) and complete the development after diapause termination in the next spring. This assumption is in line with current knowledge of the life cycle of the pest, reporting the end of post-overwintering larval development between May and July depending on local weather conditions (Potter and Held 2002;Marianelli et al. 2018;Shanovich et al. 2019). In accordance with the results of Ludwig (1928) that reported a complete inhibition of larval development during the diapause period, in our model the development of diapausing larvae is set equal to zero. Soil temperature is one of the main drivers influencing the species' diapause induction and termination (Ludwig 1932(Ludwig , 1939Fleming 1972;Potter and Held 2002). Thus, in the sub-model M1 diapause termination is triggered when the average soil temperature of the previous 15 days is above a soil temperature threshold T La to be estimated during the parameterisation process. The period of 15 days is introduced to simulate the time required for diapause termination phase and for the activation of the physiological responses related to the transition to more favourable weather conditions (Gill et al. 2017). Once diapause is terminated, individuals within the larval stage begin developing as post-overwintering larvae. The sub-model M2 simulates the development processes of post-overwintering larvae and pupae, and the flux of individuals among stages. The output of M2 is represented by the cumulative percentage of emerging adult individuals, ranging from 0 (no adult emergence occurred) to 100% (all adults are emerged).

Mathematical model
The phenology of P. japonica is simulated through a mechanistic, stage-structured demographic model based on the Kolmogorov equation (Gardiner 1985;Buffoni and Pasquali 2010). The model simulates the influence of average soil temperature on diapause termination and the nonlinear influence of soil temperature on the duration of each developmental stage (in days). The model includes a stochastic component simulating the variability of the development rates (see Section S1 of supplementary materials for details). Three developmental stages i are considered, third instar larvae ( i = 1 ), pupae ( i = 2) and adults ( i = 3 ). The physiological age x of an individual is defined as the proportion of individual development within a stage Pasquali 2007, 2010), with x = 0 representing the beginning of the development and x = 1 representing completion of the development within a stage. We introduce the condition that T La is the average soil temperature threshold, calculated in the previous 15 days, that triggers diapause termination. For the third instar larvae (i = 1) and the pupal stage (i = 2), individual physiological age varies according to a stage-specific development rate function v i (T) defined by the Brière function (Brière et al. 1999): where r i is an empirical constant, T is temperature and T i inf and T i sup are the lower and the upper developmental temperature thresholds, respectively. The development rate of the adult stage ( i = 3 ) is not defined as we consider only the emergence of the adults (i.e. the flux of individuals transferred from the pupal stage). In our model we consider T La as equal to the minimum temperature threshold allowing the development of post-overwintering larvae T La = T 1 inf . This is justified by the reasonable biological assumption that the temperature determining diapause termination and triggering larval development is equal to the minimum temperature that allows development in post-overwintering larvae. The soil temperature threshold triggering diapause termination T La and the parameters related to the development rate function of post-overwintering larvae are estimated during the parameterisation process (see Sect. 2.4.1). Parameters related to the development rate function of pupae are estimated from data available in Ludwig (1928) A graphical representation of the conceptual model and the parameters estimated in our study is presented in Fig. 1. In all the simulations, the initial conditions are set to 100 diapausing larvae at the 1st of January and to zero individuals for the other stages. All simulations end at the 31st of December of the same year. The Kolmogorov equation used for simulating pest population dynamics and phenology in sub-model M2 is presented in Gilioli et al. (2016Gilioli et al. ( , 2021, Pasquali et al. (2019Pasquali et al. ( , 2020. The full mathematical description of the model presented here can be found in Section S1 of supplementary materials.

Data on adult emergence
Time-series monitoring data between 2015 and 2019 used for model parameterisation and validation were collected by the Regional Phytosanitary Service of the Lombardy Region (Italy). A hexagonal mesh was superimposed on the area potentially infested by the species. Each hexagon has a surface of 5.42 km 2 and represents the spatial unit considered in our study (hereinafter, cell). In each cell, adult beetle activity was monitored using Tréce® traps (EPPO 2016) baited with a combination of floral attractant and synthetic pheromone lures that attract both sexes. The lures last for one season before a replacement is needed. Traps were inspected and emptied on a weekly or bi-weekly basis; then, in each cell, the adult individuals were weighed and the cell-specific population abundance was estimated considering 0.12 g as the average weight of a single individual. Within each cell, we calculated the cumulative percentage of emergence of the adult stage, as the relative cumulative adult abundance in each sampling period with respect to the total amount of adults collected in the whole flight period. As the number of traps varies among cells and among sampling years, results on adult population abundance are reported as the mean adult abundance per trap. For the purpose of our study, we have focussed on data covering the whole flight period of the adult stage (in Lombardy, approximately from May to October). Therefore, cells characterised by incomplete data were disregarded.

Environmental drivers
Data on soil temperature (depth 10 cm) at hourly temporal resolution are used as input for both M1 and M2. Although P. japonica larvae move vertically within the soil during the season (Hawley 1944), we consider soil temperature data at 10 cm depth as a suitable approximation of the overall vertical distribution of larvae and pupae during the whole period from January to May (Fleming 1972). Soil temperature datasets used in the present analysis refer to the 5th generation of European ReAnalisys (hereinafter, ERA5-land) hourly data. The ERA5-land is a reanalysis dataset which provides global data for the land component combining models with observations. This allows obtaining a globally complete and consistent dataset at a high spatial (0.1° × 0.1°) and temporal (hourly) resolution (Muñoz-Sabater 2021). The ERA5land provide soil temperature datasets at different depths. In the present study, we have extracted data from the first two soil layers (available at the average depths of 3.5 cm and 17.5 cm). Soil temperature data are then vertically interpolated using linear interpolation in order to obtain soil temperature at 10 cm soil depth. Horizontal interpolation using inverse weighting distance is used to obtain soil temperature data at the centroid of each cell considered in the present study.

Model parameters estimation
The estimation of the temperature threshold for diapause termination T La , the parameters r 1 , T 1 inf and T 1 sup of the temperature-dependent development rate function for the postoverwintering larvae (with T 1 inf = T La ) is performed via a model parameterisation process based on the solution of an optimisation problem (see Sect. 2.4.1).
For the temperature-dependent development rate function of the pupae, the parameters r 2 , T 2 inf and T 2 sup in the Brière function (Sect. 2.1.2) are estimated from data on the duration in days of the pupal stage collected by Ludwig (1928) using least square method through the function lsqcurvefit in MATLAB version R2018a. For the fitting procedure, the Fig. 1 Conceptual scheme showing the model inputs (soil temperature, temperature threshold triggering diapause termination, and development rate functions), model output (adult emergence) and the relation between the two sub-models M1 and M2. Parameters in bold related to the soil temperature threshold triggering larval diapause termination ( T La ) and to the development rate function of post-overwintering larvae ( r 1 ;T 1 inf ;T 1 sup ) were estimated through the parameterisation process. Parameters related to the development rate function of pupae ( r 2 ;T 2 inf ;T 2 sup ) were estimated from data available in Ludwig (1928) termination tolerance on the first-order optimality is set at 10 − 6 . The development of larvae during diapause is set equal to zero.

Model parameterisation and validation
From the overall dataset of time-series adult monitoring data covering the whole flight period (see Sect. 2.2.1), a first subset of 12 locations is randomly selected and used for model parameterisation (parameterisation dataset) and a second subset of 12 locations is randomly selected and used for model validation (validation dataset) (Fig. 2). In both datasets, we have ensured that cell selection includes five different sampling years (years 2015, 2016, 2017, 2018 and 2019) and four different years of first infestation (1 year, 2 years, 3 years and 4 years) defined as the period elapsed from the first appearance of the species in the cell and the observations on adult phenology used in our study. The parameterisation dataset and the validation dataset do not share observations; if the same cell is present multiple times, then the location differs in the year in which the monitoring activity was carried out. This sampling procedure allowed calibrating the model's parameters and test model's performances including the potential role of different spatial contexts, weather conditions and ages since first infestation. The random cells selected are presented in Table 1.

Model parameterisation
The model parameterisation allows estimating the soil temperature threshold T La that triggers the termination of the diapause period in M1 and the parameters r 1 , T 1 inf and T 1 sup of the development rate function of post-overwintering larvae used in M2 (parameter T 1 inf = T La ). The parameterisation process consists of minimising the mean squared error between the simulated and the observed emergence of the adult stage of P. japonica using data of the parameterisation dataset. The simulated values are computed by applying the mathematical model introduced in Sect. 2.1.2 and reported in Section S1 of supplementary materials. The minimisations are performed using the MATLAB function fmincon, based on least squares estimation, with the function tolerance set at 10 − 12 and the step tolerance set at 10 − 12 . In Table S1, the lower and the upper bounds used for the exploration of the parameters are reported. For each parameter, the range of admissible values is set slightly larger than the range of values reported in the available literature (Ludwig 1928;Fleming 1972;Régnière et al. 1981;Potter and Held 2002). Full mathematical details on the procedure used for model parameterisation is presented in Section S2 of supplementary materials. To quantify uncertainty associated with parameters' estimation of the development rate function of post-overwintering larvae, we applied to the parameterisation sample of 12 sites the Jackknife nonparametric resampling procedure (Meyer et al. 1986;Maia et al. 2000). This is an iterative procedure in which, given a starting sample of 12 elements, 12 samples of 11 elements are defined, leaving out at the n-th iteration the n-th element from the original sample. Each Jackknife sample of sites was used to estimate a development rate function. The resulting 12 development rate functions are then used to estimate the mean development rate function for the whole population and the 95% confidence interval associated with the estimated parameters. The mean development rate function obtained through the parameterisation process is then used for model validation.

Model validation
Model validation was performed using the cumulative percentage of emergence of the adult stage calculated for each location included in the validation dataset. We tested the capability of the model to simulate two fundamental aspects of the phenology of the species: (1) the overall trend of cumulative emergence distributions of the adult stage and (2) the pattern of beginning of the adult flight. The indices used for testing model's performances are different in the two cases. For the overall adult emergence, we use the mean squared error (MSE), computing the difference between simulated adult emergence ŷ i and the observed cumulative phenological curve y i (in days 2 ) The value n represents the number of sampling points available within each location considered.
For evaluating the model's capacity to simulate the beginning of adult flight, the differences (in days) between the simulated and the observed time at the 2nd, 5th and 10th percentile of adult emergence distributions are considered. Data related to the phenology of P. japonica are discrete in time as they have been collected on a weekly or bi-weekly basis. In order to allow the comparison between the simulated and the observed beginning of adult emergence, a continuous curve of the observed adult cumulative emergence should be obtained from the discrete data collected in the field. Following Curry and Feldman (1987), we use the cumulative Erlang distribution function as approximating function of the observed emergence data (see Section S3 in supplementary materials). From the cumulative Erlang function, we have then extracted the day of beginning of adult emergence calculated at the 2nd, 5th and 10th percentile and compared them with the results obtained by the model. In Figures S1-S2 we report the results of the Erlang function fitting procedure based on the time-series data estimated for both the parameterisation and the validation dataset.

Results
Adult trap catches data collected from the 24 locations included in the present study show differences in the phenology and abundance of the species (see Table S2). The time of adult emergences is reported as ordinal day, which represents the day of the year ranging between 1 and 366. Results related to the observed emergence of P. japonica adults are presented in terms of mean ± standard deviation (SD). The beginning of the first flight, calculated at the 2nd percentile of adult emergence, ranges from day 157 (June 6) to day 182 (July 1) (168 ± 6 days). The beginning of adult emergence calculated at the 5th and 10th percentile ranges from day 160 (June 9) to day 184 (July 3) (day 172 ± 6) and from day 165 (June 14) to day 187 (July 6) (day 175 ± 6), respectively. Considering only cells that were sampled for multiple years, we observe inter-annual variability in the phenology of the adults. The beginning of adult emergence calculated at the 2nd, 5th and 10th percentile of cell 4 (sampled in 2015, 2017 and 2018) is day 167 (± 4), day 171 (± 2) and day 175 (± 1), respectively. The beginning of adult emergence in cell 5 (sampled in 2016 and 2019) is day 175 (± 9), day 179 (± 6) and day 182 (± 4) for the 2nd, 5th and 10th percentile of adult emergence. In cell 12 (sampled in 2018 and 2019) the beginning of adult emergence calculated at the 2nd, 5th and 10th percentile is day 165 (± 5), day 170 (± 6) and day 174 (± 6), respectively. The time of adult peak ranges from day 178 (June 27) to day 206 (July 25), (day 192 ± 8). The total flight period of adults calculated as the difference in days between the 99th percentile and the 2nd percentile of adult emergence ranges from 40 to 70 days (56 ± 8 days). The overall adult abundance calculated as the sum of individuals captured per trap during the whole season is highly variable, ranging from 589 to 53417 individuals (15127 ± 15021 individuals). This variability may reflect the local environmental conditions and the number of years of presence of local populations after the first establishment. The model parameterisation allowed us to obtain the soil temperature threshold for diapause termination T La = 12.8 and the parameters of the development rate functions of post-overwintering larvae ( r 1 = 4.04 × 10 −5 ;T 1 inf = 12.8;T 1 sup = 29.8 ) that ensure the best fit between simulated and the observed emergence in the 12 Jackknife samples. A summary of the parameters resulting from the parameterisation process and the parameters of the development rate function of pupae fitted from the data of Ludwig (1928) is presented in Table 2. The resulting development rate functions for post-overwintering larvae (with 95% confidence interval) and pupae are presented in Fig. 3. The uncertainty related to parameters' estimate is rather low for temperatures lower than 25 °C. For temperatures between 25 and 32.6 °C, the uncertainty increases, especially in relation to the upper bound of the confidence interval. The increase in variability is due to a single Jackknife sample, in which the estimate of the maximum developmental temperature is 32.6 °C. The optimal temperature for the development of post-overwintering larvae is 25.5 °C with pupation occurring in 36 days at this constant temperature. As index for the goodness of fit for the parameterisation process, we use the averaged MSE (± SD) between simulated and observed cumulative adult emergence data considering the 12 Jackknife samples of the parameterisation dataset. The estimated MSE is 20.40 (± 2.29). The results of the parameterisation process are summarised in Table S3. The results of model validation, based on the comparisons between the observed and the simulated adult phenology, are reported in Table 3. The graphical results of the model validation procedure are shown in Fig. 4. The average MSE (± SD) between observed and simulated adult phenology considering the whole validation dataset is 26.38 (± 21.85). The model performed well in predicting the beginning of emergence of the adult stage with a mean (± SD) difference in days between simulated and observed beginning of adult emergence of 0.5 (± 4.3), 0.6 (± 3.8) and 0.5 (± 3.9) for the 2nd, 5th and 10th percentile, respectively. The model performed particularly well in predicting the observed beginning of adult emergence in 9 out of 12 locations with absolute differences in days between the simulated and the observed beginning of adult emergence lower than five days on the 2nd, 5th and 10th percentile. In three out of 12 locations, the absolute differences in days between the simulated and the observed beginning of adult emergence are lower than eight days on the 2nd, 5th and 10th percentile. The model predicted the beginning of adult emergence would occur earlier than what was observed in four (2nd and 5th percentile) and in seven (10th percentile) out of 12 locations. The difference between predicted and observed beginning of adult emergence was as high as 7 days at the 2nd percentile and 6 days at the 5th and 10th percentile, respectively. The beginning of adult emergence predicted by the model was later than observations in eight (2nd and 5th percentile) and five (10th percentile) out of 12 locations by up to seven days. Considering only cells that were sampled for multiple years, the difference in days between simulated and observed beginning of adult emergence is-2.1 (± 5.9),-2.2 (± 4.0) and-2.7 (± 2.9) for the 2nd, 5th and 10th percentile, respectively. We implemented the model on the infested area in Lombardy Region for the year 2019. The average simulated beginning of adult emergence at the 2nd, 5th and 10th percentile over the infested area is day 162 (± 6), day 166 (± 5) and day 170 (± 5).
For comparative purposes, the model is implemented using data and parameters extracted from Regniere et al. (1981). The development rate function of post-overwintering larvae is fitted through a Brière function using data on the temperaturedependent duration of the 50th percentile category of third instar larvae (see Régnière et al. (1981) for details). The temperature triggering diapause termination is set to 10 °C (Hawley 1944;Potter and Held 2002). Also in this case, we made sure that T 1 inf = T La . The obtained parameters are the following: r 1 = 2.67 × 10 −5 ;T 1 inf = T La = 10;T 1 sup = 34 . In this specific case, the average MSE (± SD) between observed and simulated adult phenology is 106.96 (± 87.1); the difference in days between simulated and observed beginning of adult emergence is-8.7 (± 4.5),-9.8 (± 3.3) and-10.4 (± 52) for the 2nd, 5th and 10th percentile, respectively. inf and T 1 sup of the development rate function of post-overwintering larvae and the temperature threshold triggering larval diapause termination T La obtained during the parameterisation process Parameters for the development rate function of the pupae r 2 , T 2 inf and T 2 sup are estimated from laboratory data extracted from Ludwig (1928) Stage Post-overwintering larvae 4.04 × 10 − 5 12.8 29.8 12.8 Pupae 8.36 × 10 − 5 11.4 39.1 -

Discussion
In this paper, we present the structure and results of a mechanistic stage-structured model simulating the role of soil temperature on the phenology of P. japonica and including the description of the post-diapause process of the species. With respect to the only other phenological model reported in literature (Regniere et al. 1981), in our work we have described the phenology of P. japonica through the Kolmogorov equation, which is a mathematical framework that allows describing the dynamical processes of stage-structured populations. The model has a stochastic component (see Section S1 of supplementary materials) that allows simulating the variability of the development rates. This is fundamental for describing the distribution of individuals in their physiological age and the distributed delay in the times of emergence in each stage. We have ensured biological realism through constraining parameters' exploration within boundaries that are based on the current knowledge about the pest. This allowed us to realistically represent the influence of soil temperature on the life-history strategies of the species. Finally, we quantified the uncertainty in the estimation of the parameters related to the development rate function of post-overwintering larvae. Higher variability in the temperature-dependent responses of third instar larvae of P. japonica is expected for temperatures higher than 25 °C (see Table S3 for details). The lower and the upper development thresholds ( T 1 inf = 12.8 and T 1 sup = 29.8 ) of post-overwintering larvae estimated through the parameterisation process are in line with data on the lower and the upper mortality thresholds (15 °C and 35 °C, respectively) reported in Ludwig (1928) and in Fleming (1972). It must be noted, however, that not all the stages tolerate these extreme temperatures. For instance, Fleming (1972) reported that larvae were able to Fig. 3 Development rate function of post-overwintering larvae estimated through the parameterisation process (red hatched area reports 95% confidence interval obtained with the Jackknife procedure) (left) and development rate function of pupae estimated using data (red asterisks) extracted from Ludwig (1928)  develop at constant temperatures only between 17.5 and 30 °C under laboratory conditions. The lower development threshold estimated in our work (12.8 °C) is higher than the value for the third instar larvae (10 °C) estimated by Régnière et al. (1981). However, it should be noted that Régnière et al. (1981) reported long developmental times at 14 °C, with pupation occurring only in the 10% of individuals after 220 days, suggesting that the substantial development of third instar larvae occurs at temperatures higher than 10 °C. In our model, the temperature threshold triggering diapause termination T La is set equal to the lower development threshold. The estimated value of 12.8 °C is higher than the temperature threshold of 10 °C proposed by some authors (Hawley 1944;Potter and Held 2002). However, in literature are also reported values similar to the one estimated from Fig. 4 Graphical results of the model validation process performed in the 12 locations of the validation dataset. In the ordinates are reported the percentage of adult emergence (from 0 to 100%). The simulated phenological curves are represented by the blue line, and the observed cumulative adult trap catches data are represented by red asterisks our data (Fleming, 1972) or even higher as for example in Ludwig (1928) that reported no development of third instar larvae at temperatures below 20 °C. The 15 days introduced as the period with the average temperature above the lower development threshold for diapause termination seems reasonable considering the delayed response of biological systems to environmental changes. In particular, diapause termination and the starting of the post-diapause development usually require irreversible changes in insect physiology that are triggered by environmental changes that act both as a stimulus for the re-activation of physiological processes and as indicators of suitable environmental conditions (Koštál 2006;Gill et al. 2017). The model performed well in predicting the overall adult phenology (average MSE 26.38 ± 21.85) and the beginning of the adult emergence calculated at the 2nd, 5th and 10th percentile.
Compared with the results obtained through the implementation of the model using data from Régnière et al. (1981), the parameterised model presented here performs well in predicting the overall adult phenology and the beginning of the adult emergence. This result justifies the work done in ensuring a thorough parameterisation of the model. Additionally, the more rigorous biomathematical approach offered by the Kolmogorov equation allows obtaining a finer representation of the phenological processes. The model also shows good capacity to interpret inter-annual temperaturedependent variability in the phenology of P. japonica. This capacity is particularly important when using the model as a tool for supporting decision-making in pest management.
The deviation of model's outputs with respect to the observations can be interpreted as the effect of environmental drivers other than temperature that might have an impact on the phenology of the species. For instance, soil moisture or food availability might play a role in larval development and overwintering period (Ludwig 1953). Additionally, the variability in the vertical distribution of larvae and their movement in soil during the pre-and post-diapause period (Fleming 1972;Villani and Nyrop 1991;Potter and Held 2002) might have an influence on individual physiology and thus on the overall phenological pattern of the species.
Further studies investigating the influence of relevant environmental drivers on the life-history strategies of P. japonica would be particularly important to improve the model parameterisation and to increase the accuracy of the model presented. This includes, for instance, setting up ad-hoc mesocosm experiments in order to observe the larval diapause termination and the post-diapause period under controlled conditions (Ludwig 1953). However, observing larval physiological changes under natural or semi-natural conditions might be complex due to difficulties linked to sampling methodology and larval vulnerability to any manipulation both in natural environment and in rearing conditions. Despite the relatively small area occupied by the pest, we observe a high variability in the beginning of the first flights (from 6 June to 8 July) and in the time of maximum abundance of adults (from 27 June to 25 July). Thus, the availability of a tool able to predict the main phenological events of P. japonica might be useful in supporting decision-makers acting at both the local level (farmers, farming cooperatives, private individuals) and at the regional level (e.g. Regional Phytosanitary Services) to the rational and efficient management of P. japonica. For instance, the model can be used at the local level (i.e. point-based simulations) for guiding the timely displacement and inspection of traps in order to get data on the emergence and the abundance of adult individuals. Similarly, the availability of information on larval diapause termination and on the beginning of their feeding activities, as well as the rate of adult emergence, which is related to their impacts and to their capacity to spread, might guide control activities and increase the efficacy of containment measures. The process-based modelling approach used in the present study assumes that the stagespecific responses to temperature estimated through the development rate functions do not depend on the thermal history of the population but they represent the influence of temperature on the physiological responses of the species. This allows implementing the model at the regional scale (i.e. area-wide simulations) to obtain maps of the phenology of the species. These maps might support decision-makers in the prioritisation of the areas of intervention, in the implementation of early-warning systems and in planning actions aimed at the wider management of P. japonica populations. Further model's expansions will aim at introducing the egg and the first and second larval instar development rate functions and the factors triggering larval dormancy period, in order to obtain a model covering the entire life cycle of the species.
Funding Open access funding provided by Università degli Studi di Brescia within the CRUI-CARE Agreement. The present work is part of the Project GESPO (CUP: E86C18002720002)-New methodologies for the integrated management of Popillia japonica (2018)(2019)(2020)(2021)(2022) funded by Regione Lombardia within the framework "Bando per il finanziamento di progetti di ricerca in campo agricolo e forestale".
Data availability Field data related to Popillia japonica adult trap catches are owned by Regione Lombardia. Other data presented in this study are available upon reasonable request.
Code availability Simulations of the phenology of Popillia japonica have been performed using a custom code.

Conflict of interest
The authors declare that they have no conflicts of interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.