Anthonomus eugenii Cano (Coleoptera: Curculionidae) is one of the most severe pests for sweet and hot varieties of pepper (Capsicum spp.; Solanaceae). The species is distributed widely, principally in Central America, but in 2013 it was detected for the first time in the Lazio region of Italy. Modelling plays a key role in reducing chemical treatments used on Capsicum spp., but reliable predictions of pest populations require adjusted tools, as well as intense knowledge of the insect's biology and its typical environment. The main goal of this work is to describe the life cycle of A. eugenii with a physiologically based model, which links the population dynamics with the environmental parameters. More specifically, this analysis focuses on the different response of the age-structured model in relation to the development rate function in input. Two methodologies to determine the best representative development rate function suitable for simulations are proposed; the first is “a priori analysis,” whereas the second is the “a posteriori analysis.” Simulations were compared with semi-field data, collected in a controlled experimental greenhouse where A. eugenii developed in varying temperature conditions. Results showed that the model used is adequate to describe A. eugenii population dynamics and highlighted how the a posteriori analysis can be essential to (i) analyze the simulation outputs, and (ii) determine the best representing development rate function, if the a priori analysis does not provide this information sufficiently clearly.
The pepper weevil Anthonomus eugenii Cano (Coleoptera: Curculionidae) represents an increasing problem for sweet and hot varieties of pepper cultivations (Capsicum spp.), not only in Central America (EPPO 2019), its land of origin, but also in Europe where this insect has been discovered recently (Speranza et al. 2014). The first Italian outbreak of the pepper weevil in 2013 (Speranza et al. 2014) was recorded in Fondi and Monte San Biagio, 2 municipalities in the province of Latina (Lazio region). The short development time of A. eugenii in favorable conditions and its destructive feeding caused considerable economic losses among pepper producers, arousing the interest of Integrated Pest Management (IPM) scientists.
Because A. eugenii may develop on both wild and cultivated Solanum and Capsicum (both Solanaceae) species, insect control strategies often are difficult to implement successfully, because nightshade-residing populations are able to re-infest cultivated pepper fields the following season (Addesso et al. 2007). Damages are inflicted primarily by the trophic activities of the larvae and adults. In addition, the oviposition behavior of adults cause abortion and fall of flowers (Patrock & Schuster 1992) and infested fruits, causing extensive losses. A complete generation usually requires 20 to 30 d (Capinera 2011), though this is strictly dependent on the optimum climatic conditions. In a laboratory, A. eugenii may reach 8 complete life cycles; however, 3 to 5 generations develop in most pepper weevil native areas (Wu et al. 2019). Indeed, it is common to find adults in the fields from Mar until Jun in the native areas of the pepper weevil, but a certain amount of adults also can be recovered throughout the yr, except for the mo of Dec and Jan (Capinera 2011).
Toapanta et al. (2005) conducted a detailed study describing the response of the pepper weevil to external environmental temperatures, which provided life tables (Harcourt 1969). An external temperature of 30 °C seems to be optimal for egg to adult development time, reproduction, fecundity, and fertility (Toapanta et al. 2005).
A mathematical description of the life cycle of pepper weevil may be of significant importance in planning control strategies to protect pepper crops. In fact, forecasts in agriculture have the potential to provide a general picture of what will happen in the following period. This feature is typical of the Decision Support System (DSS); once simulation results show a trend of an insect pest's development, experts can fix a threshold beyond which a series of consequences have to be activated in order to implement the appropriate control strategy.
The principal goal of this paper is to describe mathematically the life cycle of pepper weevil, applying the innovative approach described in the work of Rossini et al. (2019a). In particular, this study highlighted the reliability of the model applied in semi-field contexts such as greenhouses, where pepper is usually cultivated in Italy. In addition, an application protocol is proposed in order to improve future applications of physiologically based models. For this purpose, simulation results are compared with unpublished data from an experimental greenhouse located in Fondi (Latina, Italy), in the epicentre of the area that is typically infested.
Materials and Methods
POPULATION MODELLING
The mathematical details and the computational tools applied in this study are described in the works of Rossini et al. (2019a, b; 2020a). In particular, the model proposed by these authors describes an insect life cycle considering its development over time and through the life stages, and associated with the environmental temperature. The mathematical representation of the model is as follows:
where x and t indicate physiological age and time, respectively. The expressions G (x, t), β (x, t), and M (x, t) represent, in order, the development, fertility, and mortality rate functions usually estimated through laboratory experiments. The physiological age, as described in the work of Rossini et al. (2019a, 2020a), may be associated with the preimaginal stages of the life cycle of A. eugenii.
The expression G (x, t) introduces different phenological models, often represented by the development rate functions. In the case of insects, the established practice is to consider temperature as the most significant environmental variable, which drives the development of the individuals through their life stages (Mirhosseini et al. 2017; Ikemoto & Kiritani 2019). Accordingly, only the dependence of G (x, t) on time is considered. The parameters of the development rate function usually are estimated with the life tables, as described by Rossini et al. (2019b; 2020a). To distinguish the generalized form of the development rate function from its simplification, the notation changes as follows:
where T (t) is the external environmental temperature recorded in the field. In this work, the choice was made to test equation (1) for the A. eugenii life cycle with 2 particular development rate functions. The first one considered was the Briére rate function (Briére et al. 1999):
where α and m are empirical parameters, T(t) is the daily average temperature, and TL and TM are the lower and upper temperature thresholds, respectively, under and beyond which development is no longer possible for A. eugenii. The second development rate function to be considered is the Logan rate function (Logan & Hilbert 1983):
where ψ and ρ are empirical parameters, TM is the upper temperature threshold beyond which development is no longer possible, and ΔT is the interval width between the maximum peak of the function and TM .
Hence, there is a synergy between 2 apparently independent models: equation (1) manages the population dynamics, whereas equations (3) and (4) describe the insect's response to the daily average temperature. Even if equations (3) and (4) are used equally in literature, their mathematical features may provoke different responses in the simulations. In particular, the development rate functions' parameters are estimated with the same data, but the non-linear fit results may underline a preference for one expression instead of the other. In other cases, it can be difficult to highlight which is the most effective development rate function to use to represent the data. Hence, only a comparison between simulations and field data may determine these differences. This aspect will be discussed in the following section.
DATA ANALYSIS
Once the non-linear fit results have provided the functions' parameters, the choice of the best representative development rate function between equations (3) and (4) will be determined with a χ2-test and by analyzing the adjusted R2-value. This procedure, together with computational aspects, was described by Rossini et al. (2019a). To estimate the parameters of the expressions (3) and (4) the ROOT's “Minuit” tool (Brun & Rademakers 2007) was included in an ad hoc C++ program, as reported by Rossini et al. (2020a). As mentioned previously, it may happen that the and values resulting from the non-linear regression with equations (3) and (4) are similar. When this condition is verified, an additional operation is required to assess which is the best representative development rate function for the specific application. In a preliminary analysis, if the non-linear fits results have minor or otherwise negligible differences, the user could make a direct choice of which development rate function to use. This statement is supported by the empirical formulation of both Briére and Logan formulas.
On the other hand, an alternative is to run simulations considering both the Briére and Logan expressions as development rate functions as input into equation (1). Simulation outputs then will be compared with field data, and the most effective development rate function to use in further applications will be determined a posteriori. The comparison between simulation outputs and field data will be conducted using the methodology described in Rossini et al. (2019a; 2020a, b). The authors reported that the distance between the simulated and the actual experimental populations would be evaluated, again considering the function. As an indicator of the distance, a lower χ2-value means a more accurate overlap between simulations and field data. Accordingly, if the best development rate function to use as input in equation (1) cannot be selected using the criteria above, it can be helpful to compare simulation outputs with field data. Thereafter, the development rate function that provides a lower χ2-value will be applied in further simulations.
EXPERIMENTAL DESIGN FOR VALIDATION
The model validation was conducted during the 2014 growing season in an experimental greenhouse located in Fondi (Latina, Italy) at the company GENISTA S.r.l., which specializes in genetic improvement. The greenhouse had a length of 40 m, a width of 10 m, a minimum height of 1.95 m, and a maximum height of 3.5 m at the center of the roof. Pepper plants were transplanted into the greenhouse in Aug, which was opened to promote A. eugenii infestations.
Once the presence of pepper weevils was ascertained inside the experimental field by visual observations, the greenhouse was isolated with an anti-aphid net (Agritech s.r.l., Eboli, Salerno, Italy) to avoid the further entrance or exit of A. eugenii individuals. At this point, the adults enclosed in the greenhouse mated and laid eggs on pepper plants, producing the generation that was monitored to validate and analyze equation (1).
The monitoring activity focused on pepper weevil adults. For this purpose, 5 yellow chromotropic sticky traps (Serbios s.r.l., Badia Pole-sine, Rovigo, Italy) were deployed along the diagonal axis of the field, and checked each wk from Oct to Dec. The experiment was concluded at the end of Dec when the pepper fruit was harvested.
The environmental temperature in the greenhouse was recorded from the date the pepper plants were transplanted until the interruption of the experiment. More specifically, a meteorological station placed in the greenhouse provided the daily average temperature based on 24 measurements (1 every h).
The error associated with the experimental population (dots in Figs. 1 & 2) is related to sampling time. Normally, this value is fixed at half the temporal range between the 2 field surveys, which is ± 3.5 d.
Results
The life tables provided by Toapanta et al. (2005) facilitated the estimation of the development rate function parameters. The best-fit parameters for equations (3) and (4) are listed in Table 1, together with χ2 and R2 values. Figures 3 and 4 provide a graphical representation of the best-fit functions. Additional information (correlation matrix, variances and covariances matrix, and a plot of the adjustment curves with 95% confidence intervals) are reported in the supplementary materials.
The criteria for selecting the best development rate function to use for simulations, as described in the “Data Analysis” section, requires a cross check between the χ2-test and the R2-value, if the and values are not too similar between the 2 functions. The values listed in Table 1 underscore that this criteria cannot be applied, since a χ2-test indicates a P < 0.01 for both equations (3) and (4). In addition, the respective coefficients of determination are similar. On the basis of the χ2-test, there is a preference for the Briére equation, but considering the R2, the situation is overturned. These first results highlight that both Logan and Briére equations represent the life tables from Toapanta et al. (2005) effectively; nevertheless, it remains difficult to choose 1 of the 2 equations. From a biological point of view, it could be more interesting to select the Briére equation instead of the one by Logan. The reason for this is that the information about the lower temperature threshold for the development of the species, T_L, is not indicated by the Logan equation. On the other hand, the aim of this work is to provide an adjusted tool to use as a decision support system; hence, the best representative development rate function will be determined a posteriori, comparing the simulations with field data.
Fig. 1.
Simulation output from the model (1) evaluating the daily average temperature with the Briére development rate function.
An additional and important piece of biological information provided by this first part of the analysis concerns the optimal thermal requirements for the development of A. eugenii. In particular, an optimal temperature of (31 ± 1) °C has been calculated for the Briére equation, whereas (29 ± 3) °C has been calculated for the Logan equation. Considering the errors associated with the optimal temperatures provided by the development rate function equations, it is possible to determine that these results are, in fact, in accordance with each other.
Fig. 2.
Simulation output from the model (1) evaluating the daily average temperature with the Logan development rate function.
Table 1.
The parameter values for Briére and Logan rate functions reported with their standard errors, calculated with the methodology reported by Rossini et al. (2019a, b; 2020a), using data from Toapanta et al. (2005). Additional information includes the χ2-value, the one-tailed probability of the distribution, R2-value, and the number of degrees of freedom (NDF).
Simulation results are reported in Figures 1 and 3. In particular, Figure 1 reports the simulation output compared with field data considering the Briére equation as input. In this case, the projections indicate a maximum peak of the population on 22 Nov, whereas field data confirmed that the maximum occurred on 23 Nov. The same situation was reported for the Logan equation, as shown in Figure 4.
The principal differences between the 2 simulations are concentrated in the right-hand side of the plot, on the decreasing tails of the graph. Table 2 reports the χ2 and R2 values and: whereas the first value indicates the distance between simulations and field data, the second value is helpful for understanding the extent to which the simulations represent the actual greenhouse population. On the basis of the values in Table 2, it is possible to assess that simulations are more reliable if, in equation (1), the response of A. eugenii to the daily average temperature is evaluated by the Logan equation.
Discussion
The results obtained highlight that the Logan equation (inserted into equation [1]) is more appropriate to describe A. eugenii populations. This is not true for all insect species, and it depends on several aspects. For example, the non-linear fit results with equations (3) and (4) are strongly dependent on the number of experimental points in the life tables. Accordingly, enriching the life tables with additional points can increase the accuracy of the first part of the analysis. Besides demonstrating and analyzing these results, this work aims to provide an operating protocol to apply physiologically based models, which requires a combination of a population dynamics model and a phenological model. Specifically for this work, population dynamics is managed by equation (1), while the phenology (in this case, the relationship between insects and environmental temperature) is managed by the development rate function expressions.
Fig. 3.
Briére development rate function compared with life tables point from Toapanta et al. (2005).
An additional result, as demonstrated above, is the utility of the a posteriori analysis. It has, in fact, been shown that when it is impossible to decide prior to the simulation the most appropriate development rate function to use as input in the population dynamics model, the a posteriori analysis is the only viable means of reaching a decision in this respect.
As already anticipated, both Briére and Logan equations are formulated on an empirical basis; there are no reasons to prefer one instead of the other because they represent the life tables in a similar way. On the other hand, the a posteriori analysis showed that even if they represent the life tables in a similar way, the simulation outputs can report more discernible differences.
The considerations resulting from this study are helpful for model scientists who aim to provide tools for direct applications in integrated pest management. The use of mathematical models as a decision support system is growing, and its application has been found in several case studies (Rupnik et al. 2019), most of which are sponsored by government programs (Carberry et al. 2002). For instance, there are models which consider not only the potential abundance of pests in the field, but also include population density in relation to certain economic thresholds (Tang & Cheke 2008). Despite the existence of various proposals for the mathematical representation of population dynamics (Nance et al. 2018; Ainseba et al. 2011; Rossini et al. 2019b, 2020b), the reliability of the simulations is strongly dependent on the function which relates the species with its environmental parameters. Accordingly, it is not the best method to choose a priori a specific development rate function. Instead, the development rate functions should be evaluated twice: (i) a priori to exclude wide differences in representing life table data, and (ii) a posteriori to evaluate potential differences between different functions in input.
Table 2.
Numerical evaluation of the differences between the simulation outputs and the field data. In this case, χ2 indicates the distance between simulations and field data, whereas R2 indicates how much the simulation represents field data.
Acknowledgments
The authors thank the reviewers for their comments and suggestions, which have been helpful for the improvement of this manuscript. The authors gratefully acknowledge Carlo Vagnozzi and Jessica Santianni from GENISTA S.r.l. for their kind concession of the experimental greenhouse, and Simone Pesolillo and Emanuele Eusepi from Meteotec S.r.l. for the climatic data acquisition and management. The research was carried out in the framework of the Ministry for Education, University and Research (MIUR) initiative “Department of Excellence” (Law 232/2016). “Pilot project on the formation of pepper weevil-free areas and calibration of a mathematical model to predict the severity of pepper infestation in the Fondi area.” Funded by Reg. 1698/2005 PSR 2007/2013 of Lazio Region. Misure 124, Cooperation for the development of new products, processes and technologies in the agricultural, food and forestry sectors. Public Notice DGR n. 76/2014. Act granting aid N. 28/124/10 of 18/12/2014, n. 8475921005.
