A Temperature-Dependent Model for Tritrophic Interactions Involving Tea Plants, Tea Green Leafhoppers and Natural Enemies

Simple Summary Tritrophic interactions have achieved much attention in research on ecology. The tea plant Camellia sinensis (L.) O. Kuntze is a major economic crop in Asian countries, especially in China. Tea plants suffer infestations from herbivory attack during their lifetime. The tea green leafhopper (Empoasca onukii Matsuda) is a major pest for tea plants. The parasitic or predatory natural enemies of tea green leafhoppers (TGLs) can feed on their eggs, nymphs, or adults. However, a detailed mathematical model for tea plants–TGLs–natural enemies is still lacking. In the current work, we established a novel model based on laboratory measurements or field observations with temperature-dependent effects on tritrophic interactions for tea ecosystems. As expected, cyclic behaviors are identified. Stochastic simulations further showed two TGL outbreaks, the timing of which is consistent with field observations. Effective accumulated temperature (EAT) is possibly an important predictor of TGL outbreak. Applying slow-releasing semiochemicals as either repellents or attractants may be highly efficacious for pest biocontrol. An optimal treatment time of semiochemicals can also be determined. Our detailed model identifies key features of tritrophic interactions involving tea plants and can be extended to other ecosystems. Abstract The tea green leaf hopper, Empoasca onukii Matsuda, is a severe pest of tea plants. Volatile emissions from tea shoots infested by the tea green leafhopper may directly repel insect feeding or attract natural enemies. Many studies have been conducted on various aspects of the tritrophic relationship involving tea plants, tea green leafhoppers and natural enemies. However, mathematic models which could explain the dynamic mechanisms of this tritrophic interaction are still lacking. In the current work, we constructed a realistic and stochastic model with temperature-dependent features to characterize the tritrophic interactions in the tea agroecosystem. Model outputs showed that two leafhopper outbreaks occur in a year, with their features being consistent with field observations. Simulations showed that daily average effective accumulated temperature (EAT) might be an important metric for outbreak prediction. We also showed that application of slow-releasing semiochemicals, as either repellents or attractants, may be highly efficacious for pest biocontrol and can significantly increase tea yields. Furthermore, the start date of applying semiochemicals can be optimized to effectively increase tea yields. The current model qualitatively characterizes key features of the tritrophic interactions and provides critical insight into pest control in tea ecosystems.


Introduction
The tea plant Camellia sinensis (L.) O. Kuntze is a major economic crop in Asian countries, and tea is also a non-alcoholic drink worldwide [1]. As with other crops, tea plants suffer infestations from herbivores during their lifetime. Specifically, the tea green leafhopper (TGL), Empoasca onukii Matsuda (Hemiptera: Cicadellidae), is probably the most serious constraint to tea plant cultivation [2]. Nymphs or adults of TGLs attack the tea plants by sucking the sap, and their outbreaks can significantly reduce tea yields and quality [3].
Plants have evolved diverse mechanisms against herbivore attack [4]. Herbivory attack induces the release of plant volatiles, termed as herbivore-induced plant volatiles (HIPVs) [4]. These HIPVs are involved in tea plant defense, either directly by repelling herbivores, or indirectly by attracting natural enemies of herbivores [5,6]. For example, two parasitic wasps, Stethynium empoascae Subba Rao and Schizophragma parvula Ogloblin can parasitize the eggs of tea green leafhoppers, with a parasitism rate of 40 to 65% [2,7,8]. Recent reports have shown that HIPVs induced by tea green leafhopper infestation in tea plants can attract these parasitic wasps, which serve as natural enemies of tea green leafhoppers [9,10]. This volatile-mediated tritrophic interaction between tea plants, tea green leafhoppers and parasitoid wasps has received much more attention in integrated pest management [11,12].
In the past few decades, interest in the design of mathematical models based on population dynamics to aid pest management has increased [13][14][15][16][17]. Most studies which focused on predator-prey interactions often assume that the environments are homogeneous [18,19]. The natural system, however, is inevitably influenced by various environmental factors (e.g., temperature). Temperature is a determinant factor and has significant impact on processes related to ecological community, particularly predator-prey or consumer-resource interactions [20,21]. For example, the growth rate of crops, the oviposition, and the lifespan of pests can be altered by temperature [22,23]. For predator-prey systems, variations in temperature usually have a nonlinear effect on the short-term interaction strength by regulating the foraging and attack rates of predators [24]. Furthermore, the thermal effect on prey traits is also relevant to the predator-prey interactions [25]. Many species, including insects, undergo ontogenetic shifts, either in the form of habitats or life cycles with distinct developmental stages [20]. This is critical, since the consumption of some prey species is dependent on the prey body size or its developmental stage [26]. Predator and prey traits may sometimes have asymmetric thermal responses; for instance, the optimal temperatures for foraging and development are different, or predator and prey traits respond differentially across a thermal gradient [20]. Under these circumstances, dynamic behaviors in tritrophic or interlocked predator-prey systems may become complicated owing to temperature-dependent effects.
Additionally, there is considerable evidence that noise or randomness can play an important role in ecological systems [27]. These noisy fluctuations can arise from either intrinsic sources, which are associated with inherent stochasticity in ecological models (e.g., randomness in prey selection), or from extrinsic sources, which correlate with environmental factors (e.g., temperature) [28]. As a result, deterministic models are imprecise for characterizing the dynamic evolution of tritrophic interactions in tea agroecosystems, and stochastic models should be used. However, to the best of our knowledge, a temperaturedependent stochastic model which characterizes the tritrophic interactions among tea plants, tea green leafhoppers, and natural enemies is still lacking. To facilitate better understanding of tritrophic interactions in tea ecosystems, we constructed a detailed and realistic stochastic model based on ordinary differential equations (ODEs) with variations in daily temperature and thermal effects on tea growth and developmental stages of tea green leafhoppers. Stochastic simulations identified two peaks or outbreaks of tea green leafhoppers in a year, with markable variations in the peaking time. The role of parasitoids or predators in the tritrophic interactions was evaluated. The results indicated that predators may be highly efficacious for leafhopper control compared with parasitoids. The relationship between daily average effective accumulated temperature (EAT) and tea green leafhopper outbreaks was also established. An optimum treatment using semiochemicals was also proposed to increase tea yields. Collectively, our model characterized the dynamic thermal responses of tritrophic relationship involving tea plants, tea green leafhoppers, and natural enemies to provide important insight into the development of optimum strategies for tea pest control.

Model Construction
We first constructed a tea plant-tea green leafhopper-natural enemy model, and the output was used to propose optimum biocontrol strategies. In model (1), C represents the tea crop population, whereas P and W denote adult tea green leafhoppers and parasitic wasps, respectively. E denotes leafhopper eggs and N represents the nymph population. Previously, Holling proposed Holling-type functional responses in the predator-prey system [29]. Recently, several authors have found that the Holling type II response can make the standard predator-prey model more realistic [30,31]. Therefore, we assumed that the consumption of adult leafhoppers or nymphs follows Holling type II functional responses [29]. r is the intrinsic growth rate of tea plants and K is the environmental capacity. Since the adults and nymphs of tea green leafhoppers can pierce and suck the sap of tender tea shoots [32], P and N can both feed on tea crops. a 1 denotes the predation rate of pests. a 2 indicates the parasitic rate. h is the harvesting or death rate of tea crops. c 1 represents the conversion efficiency from adult leafhoppers to eggs by feeding on tea crops. Similarly, c 2 is the conversion efficiency of parasitic wasps from the egg population. b 1 and b 2 are the half saturation rates for the feeding process. µ represents the natural death rate of pests. δ describes the dispersal rate for parasitic wasps (flying away from the sites of infestation). The parameters were generally adopted from the work by Mandal et al. with a few modifications to fit our tea plant-tea green leafhopper-natural enemy system or observations [31]. The feeding damage on the tea crop by leafhoppers induces the synthesis and emission of VOCs which can attract natural enemies ( Figure 1A) [5,6]. λ characterizes the attraction rate of parasitic wasps induced by feeding damage. The developmental stages of leafhoppers include eggs, nymphs, and adults [33]. ω 1 is the (average) conversion rate from eggs to nymphs, and ω 2 is the (average) conversion rate from nymphs to adult leafhoppers. The following ODEs were developed based on the above assumptions. For parameter values and initial conditions, please refer to Table S1.  TGLs induces the release of volatiles, which can attract the wasps. (B) Limit cycles in 3D plot with increasing growth rate of tea crops (r). In the simulation, r is increased from 4 to 12, with an increment of 2. (C) Limit cycles in 3D plot with increasing parasitic rate (a2). a2 is increased from 0.4 to 2, by 0.4. (D) Local sensitivity analysis for period (blue), amplitude (red), and integrated tea crop levels (green, integral, or tea yields).

Stochastic Tea Pest-Natural Enemies Model with Thermal Effect
The deterministic model only describes average population dynamics and may be imprecise if stochasticity is considered. We further constructed a realistic and stochastic model which incorporates the effect of daily temperature on tea crop growth, developmental stages, and lifespan of tea green leafhoppers. The effect of temperature on plant growth rate was formulated as previously described [34]. Briefly, the coefficient for temperature-driven effect (TF) can be expressed as: where T0 denotes the optimal temperature on tea growth. Following the work by Chen et al. [23], T0 was set to be 23 °C, and β was estimated to be 20 °C. The temperature-dependent tea growth rate can then be expressed as: The tea growth rate r was replaced by r(T) and modified in stochastic simulations. The egg, nymph, and pre-oviposition stages during tea green leafhopper development are also temperature-dependent [35]. The relationship between temperature and stage durations can be described using the Arrhenius equation bT a e − ⋅ [20], where T is the temperature. The parameter fitting results are listed in Table S2 and Figure S1. Note that with increasing growth rate of tea crops (r). In the simulation, r is increased from 4 to 12, with an increment of 2. (C) Limit cycles in 3D plot with increasing parasitic rate (a 2 ). a 2 is increased from 0.4 to 2, by 0.4. (D) Local sensitivity analysis for period (blue), amplitude (red), and integrated tea crop levels (green, integral, or tea yields).

Stochastic Tea Pest-Natural Enemies Model with Thermal Effect
The deterministic model only describes average population dynamics and may be imprecise if stochasticity is considered. We further constructed a realistic and stochastic model which incorporates the effect of daily temperature on tea crop growth, developmental stages, and lifespan of tea green leafhoppers. The effect of temperature on plant growth rate was formulated as previously described [34]. Briefly, the coefficient for temperaturedriven effect (TF) can be expressed as: where T 0 denotes the optimal temperature on tea growth. Following the work by Chen et al. [23], T 0 was set to be 23 • C, and β was estimated to be 20 • C. The temperaturedependent tea growth rate can then be expressed as: The tea growth rate r was replaced by r(T) and modified in stochastic simulations. The egg, nymph, and pre-oviposition stages during tea green leafhopper development are also temperature-dependent [35]. The relationship between temperature and stage durations can be described using the Arrhenius equation a · e −bT [20], where T is the temperature. The parameter fitting results are listed in Table S2 and Figure S1. Note that in this case, the logistic function was not statistically significant, as all fitted parameters had confidence intervals that contained zeros. The effect of temperature on mean lifespan of eggs, nymphs, and adults was also fitted using the Arrhenius equations. The fitting results are shown in Table S2.
According to the law of exponential decay, the kinetic parameters for degradation/ conversion rates σ and the mean lifetime/duration of a stage τ can be interchanged using a reciprocal relation.
Therefore, µ, ω 1 , or ω 2 can be regarded as degradation/conversion rates and are functions of temperature T. The parameters µ, ω 1 , or ω 2 were replaced by µ(T), ω 1 (T), or ω 2 (T) with 10% variations. We did not differentiate the five instars during the nymph stage. The temporal duration of nymph and pre-oviposition stages were combined to calculate ω 2 (T), because leafhoppers at nymph and pre-oviposition stage can both feed on tea crops but are immature (unable to lay eggs). The threshold for the different developmental stages was also considered. The threshold temperatures determined for the egg and nymph (preoviposition) stages are approximately 11.06 • C and 7.86 • C, respectively [35]. Therefore, we obtained the following expressions: T θ 1 is the developmental threshold temperature for eggs. T θ 1 was set to be 11 • C in our model. T θ 2 is the threshold temperature for nymphs. T θ 2 was 8 • C in our model. The final equations for stochastic simulations are shown as follows: The temperature (T) varies over a year and, therefore, we obtained the average daily temperature at Anqing City in 2020 from the National Oceanic and Atmospheric Administration website ( Figure S3, http://www.noaa.gov/, accessed on 5 May 2022).

Local Sensitivity Coefficient
Suppose a dynamic system x' = F(x, θ), where x and θ are state and parameter vector, respectively. θ = (θ 1 , θ 2 , . . . , θ i , . . . , θ n ), where n is the number of model parameters. For a certain feature 'I' (period, amplitude, or integral), its local sensitivity coefficient S with respect to changes in parameter θ i is defined as: When 'I' stands for integrated response or integral, it is defined as: where T end is the end of a year, and A denotes the level of certain population (tea crop, eggs, adult leafhoppers, or nymphs). The integrated response of tea crops is regarded as tea yields. In numerical simulations, the calculation of local sensitivities is obtained by increasing the parameter θ i by 1% and keeping all the other parameters fixed; we then get the perturbed response I p for feature 'I'. I default is response for feature 'I' using default parameters (Table S1). Thus, in our model, ∂I/I = (I p − I de f ault )/I de f ault , and ∂θ i /θ i = 1%.

Stochastic Simulation
Stochastic simulation based on ODEs was implemented using a Poisson τ-leap method [36]. Briefly, suppose there are N populations in the system [S 1 , S 2 , . . . , S i , . . . , S N ] with molecular number x = X i (t) for species S i at time t. The reaction channels (list of reactions) are represented as [R 1 , R 2 , . . . , R M ], and M is the number of reactions. For the reaction channel R j , a j (x)dt denotes the probability that the reaction R j will occur in the time interval [t, t + dt). A state change vector ν j is defined to characterize reaction channel R j . ν ij is the element of ν j and represents the change in molecular numbers of population S i due to the reaction R j .
In the Poisson τ-leap method, we assume that there are many reactions firing in a relatively large time interval [t, t + τ) [36]. The number of reactions for channel R j firing at the interval [t, t + τ) obeys Poisson distribution with mean a j (x)τ. The system is updated by: 'P' denotes the Poisson distribution. The basic framework for implementation of this stochastic algorithm based on ordinary differential equations is provided in Supplementary Materials.

Model Simulation
The ordinary differential equations were integrated with the ode23s solver. Numerical and stochastic simulations were performed in MATLAB (MathWorks, Natick, MA, USA, R2018b).

The Tea Plant-TGL-Wasp System Exhibits Sustained Oscillation
The tritrophic interactions among tea plants, TGLs, and parasitic wasps were formulated using model 1. After infestation, the feeding damage by TGLs induces the emission of volatiles from the tea crop, and the VOCs may attract natural enemies (dominant parasitoids S. empoascae or inferior species S. parvula) ( Figure 1A). Using the default parameters, the system exhibited sustained oscillations, as evidenced by limit cycles ( Figure 1B,C and Table S1). We noted the existence of interlocked positive (tea-parasitic wasps-leafhoppers) and negative (tea-leafhoppers and leafhoppers-parasitic wasps) feedback loops in the tritrophic system, which may potentially lead to limit cycle oscillations [37]. Indeed, a supercritical Hopf bifurcation can be identified ( Figure S2). The attraction of natural enemies to the sites of herbivory attack by volatiles creates a positive feedback loop in the tritrophic interaction system (Figure 1A), and the positive feedback may function to stabilize the limit cycle oscillations [37]. Increasing the growth rate of tea plants enlarged the limit cycles ( Figure 1B), whereas increasing the parasitic rate (a 2 ) slightly shrank the size of the limit cycles ( Figure 1C). Local sensitivity analysis showed that parameters for tea growth (r and K) negatively affected the period, whereas enhanced insect-feeding processes (a 1 and b 1 ) tended to elongate the period ( Figure 1D). The process related to tea growth also increased the amplitude and temporal integral (tea yields) of tea plants ( Figure 1D). TGL feeding indeed decreased the tea yields ( Figure 1D). We then extended the parametric regions and evaluated the systemic responses. Undamped oscillations were allowed in broad regions for a 1 and λ (Figure 2A,B). Interactions between a 1 and λ were also identified.
For relatively large λ values, the effect of a 1 on oscillatory periods and amplitudes was weakened (Figure 2A,B). There was also an interaction between parameters r and a 2 , as we observed that the effect of r on oscillatory patterns vanished when the value of a 2 was large ( Figure 2D,E). When sustained oscillations were dampened in some regions (lower a 1 /higher λ, or higher a 2 and r), the tea crop integral (tea yields) was strongly increased ( Figure 2C,F). These results suggest that the tritrophic relationship among tea plants, TGLs, and parasitic wasps allows the system to perform cyclic behavior, and that regulation of key parameters can strongly affect tea yields.
that parameters for tea growth (r and K) negatively affected the period, whereas enhanced insect-feeding processes (a1 and b1) tended to elongate the period ( Figure 1D). The process related to tea growth also increased the amplitude and temporal integral (tea yields) of tea plants ( Figure 1D). TGL feeding indeed decreased the tea yields ( Figure 1D). We then extended the parametric regions and evaluated the systemic responses. Undamped oscillations were allowed in broad regions for a1 and λ (Figure 2A,B). Interactions between a1 and λ were also identified. For relatively large λ values, the effect of a1 on oscillatory periods and amplitudes was weakened (Figure 2A,B). There was also an interaction between parameters r and a2, as we observed that the effect of r on oscillatory patterns vanished when the value of a2 was large ( Figure 2D,E). When sustained oscillations were dampened in some regions (lower a1/higher λ, or higher a2 and r), the tea crop integral (tea yields) was strongly increased ( Figure 2C,F). These results suggest that the tritrophic relationship among tea plants, TGLs, and parasitic wasps allows the system to perform cyclic behavior, and that regulation of key parameters can strongly affect tea yields.

Stochastic Simulation Identifies Properties of TGL Outbreaks
We next performed stochastic simulations with temperature-dependent parameters. The daily temperature in 2020 was used as an example ( Figure S3). In contrast to the deterministic simulations using ODEs, stochastic trajectories with temperaturedependent regulation indeed showed highly variable dynamics ( Figure 3A,B). There were generally two outbreaks of tea green leafhoppers in one year, together with two peaks in the number of eggs ( Figure 3A,B). Estimation of the peaks of pest outbreaks showed two skewed distributions ( Figure 3C). The median date of the first outbreak was 24 May with 95% confidence interval (CI) [17 May,31 May], whereas the median for the second outbreak was 30 September [1 September, 29 October] ( Figure 3C). In our recent field experiments, two leafhopper outbreaks were identified ( Figure 3D). The estimated time for outbreaks in the model was consistent with our field observations ( Figure 3C,D) and previous studies [38,39]. In addition, the number of eggs also had a trough, with a median on 6 July with a 95% CI [30 June, 12 July] ( Figure 3E), which is compatible with observations [40]. These simulations suggest that our stochastic model is generally consistent with field observations and shows significant variations in temporal trajectories.

Stochastic Simulation Identifies Properties of TGL Outbreaks
We next performed stochastic simulations with temperature-dependent parameters. The daily temperature in 2020 was used as an example ( Figure S3). In contrast to the deterministic simulations using ODEs, stochastic trajectories with temperature-dependent regulation indeed showed highly variable dynamics ( Figure 3A,B). There were generally two outbreaks of tea green leafhoppers in one year, together with two peaks in the number of eggs ( Figure 3A,B). Estimation of the peaks of pest outbreaks showed two skewed distributions ( Figure 3C). The median date of the first outbreak was 24 May with 95% confidence interval (CI) [17 May,31 May], whereas the median for the second outbreak was 30 September [1 September, 29 October] ( Figure 3C). In our recent field experiments, two leafhopper outbreaks were identified ( Figure 3D). The estimated time for outbreaks in the model was consistent with our field observations ( Figure 3C,D) and previous studies [38,39]. In addition, the number of eggs also had a trough, with a median on 6 July with a 95% CI [30 June, 12 July] ( Figure 3E), which is compatible with observations [40]. These simulations suggest that our stochastic model is generally consistent with field observations and shows significant variations in temporal trajectories.

Efficiency of Predators Is Higher Compared with Parasitoids in Pest Control
We next compared the efficiency of predation and parasitism for tea pest control. Equation (6) was slightly modified to incorporate the predators (Supplementary materials: Stochastic model for tea plant-TGLs-predator relationship). All parameters remained the same in both the parasitoid and predator models. When only predators were present, the tea crop integral (tea yields) was significantly higher ( Figure 4A). Predators were more effective in reducing the leafhopper population ( Figure 4B,C). As expected, the egg population recorded in the parasitoid-only model, was remarkably lower than that in the predator-only model ( Figure 4D). Since recruited predators feed on adults or nymphs of tea green leafhoppers, the predators can provide direct or instant protection against herbivory attack compared with parasitic wasps, whose protective effect is indirect and delayed. These results suggest that predators may provide better protection against TGL infestation compared with parasitic wasps, and are consistent with field observations [41].

Efficiency of Predators Is Higher Compared with Parasitoids in Pest Control
We next compared the efficiency of predation and parasitism for tea pest control. Equation (6) was slightly modified to incorporate the predators (Supplementary materials: Stochastic model for tea plant-TGLs-predator relationship). All parameters remained the same in both the parasitoid and predator models. When only predators were present, the tea crop integral (tea yields) was significantly higher ( Figure 4A). Predators were more effective in reducing the leafhopper population ( Figure 4B,C). As expected, the egg population recorded in the parasitoid-only model, was remarkably lower than that in the predator-only model ( Figure 4D). Since recruited predators feed on adults or nymphs of tea green leafhoppers, the predators can provide direct or instant protection against herbivory attack compared with parasitic wasps, whose protective effect is indirect and delayed. These results suggest that predators may provide better protection against TGL infestation compared with parasitic wasps, and are consistent with field observations [41].

Daily Average Effective Accumulated Temperature Influences the Pest Outbreaks
We next evaluated the effect of temperature on the pest outbreaks. Effective accumulated temperature (°C•day) is defined as where Ti is the (average) daily temperature of the ith day, n is the duration of record, and B is the developmental threshold temperature. B is set to be 7.8 °C [35]. The daily average EAT is obtained by The EATs were recorded for various durations (n = 10, 50, 80, 100, and 120 days). Then, 5000 stochastic simulations were performed. When pest outbreaks or peaks were identified, the daily average EAT on that day was also calculated. The number of incidences of a pest outbreak or peak on a certain day was normalized to obtain the frequency. We found that occurrence of pest outbreaks was only restricted within narrow ranges of the daily average EAT (roughly 8~23 °C, Figure 5). Daily average EATs outside this range failed to identify a pest outbreak. When the duration of record was longer, the distribution became more concentrated (n ≥ 80, Figure 5). These simulations demonstrate that the daily average EAT strongly influences the onset of pest outbreaks.

Daily Average Effective Accumulated Temperature Influences the Pest Outbreaks
We next evaluated the effect of temperature on the pest outbreaks. Effective accumulated temperature ( • C·day) is defined as where T i is the (average) daily temperature of the ith day, n is the duration of record, and B is the developmental threshold temperature. B is set to be 7.8 • C [35]. The daily average EAT is obtained by 1 The EATs were recorded for various durations (n = 10, 50, 80, 100, and 120 days). Then, 5000 stochastic simulations were performed. When pest outbreaks or peaks were identified, the daily average EAT on that day was also calculated. The number of incidences of a pest outbreak or peak on a certain day was normalized to obtain the frequency. We found that occurrence of pest outbreaks was only restricted within narrow ranges of the daily average EAT (roughly 8~23 • C, Figure 5). Daily average EATs outside this range failed to identify a pest outbreak. When the duration of record was longer, the distribution became more concentrated (n ≥ 80, Figure 5). These simulations demonstrate that the daily average EAT strongly influences the onset of pest outbreaks.

Slowly Releasing Semiochemicals Increases Tea Yields
Pheromones and other semiochemicals are widely used to manage agricultural or forest pests [42]. Synthetic semiochemicals can either be used as repellents or attractants [42]. To simulate the effect of semiochemical application with different releasing rates or durations, we simulated the effect of controlled release. The amount or dosage of semiochemicals was fixed (720 in model 1), whereas the releasing duration was varied ( Figure 6A, right). The effect of attractants was evaluated. Attractants can recruit the natural enemies of TGLs [42], analogous to an increase in λ values. We modified λ = 0.3 (1 + vmax*SCs/(SCs + Km)); SCs represent the current level of semiochemicals. The start time of semiochemical usage was either on 19 February (50th day) or on 29 April (120th day) in the simulations for case studies ( Figure 6A). A variable duration of semiochemical release was generated following a normal distribution D~N(60, 20), where D represents duration in days ( Figure 6A, right). Since the dosage of semiochemicals was fixed (D•SCs = 720), the inner areas of each square were all the same (green, yellow to red, Figure 6A, right). Variations in the durations of release resulted in changes in the efficiency of attracting wasps or parameter λ ( Figure 6A, left). Stochastic simulations were then performed, and the relation between the tea yields and duration of semiochemical release was determined. Significant correlation was identified irrespective of the starting time ( Figure 6B,C, top). Increasing the duration of releases raised the tea yields in a year ( Figure  6B,C, bottom panels, p < 0.0001). Qualitatively similar results were obtained if semiochemicals were used as repellents ( Figure S4).
We further investigated whether there existed an optimum time for the application of semiochemicals. The starting time was varied from the 30th day to the 130th day. The averaged tea yields after using semiochemicals were then determined (Supplementary Materials: Average daily tea yields after semiochemical application). We found that the averaged tea yields first increased and then decreased when the starting time was changed from the 30th day to the 130th day ( Figure 7A-H). There was a local maximum around the 90th day ( Figure 7I). These results suggest that semiochemicals usage can be optimized to obtain higher average tea yields.
where T i is the (average) daily temperature of the ith day, and B is the baseline or developmental threshold temperature. The duration of record 'n' is varied from 10, 50, 80, 100, to 120 days.

Slowly Releasing Semiochemicals Increases Tea Yields
Pheromones and other semiochemicals are widely used to manage agricultural or forest pests [42]. Synthetic semiochemicals can either be used as repellents or attractants [42]. To simulate the effect of semiochemical application with different releasing rates or durations, we simulated the effect of controlled release. The amount or dosage of semiochemicals was fixed (720 in model 1), whereas the releasing duration was varied ( Figure 6A, right). The effect of attractants was evaluated. Attractants can recruit the natural enemies of TGLs [42], analogous to an increase in λ values. We modified λ = 0.3 (1 + v max *SCs/(SCs + K m )); SCs represent the current level of semiochemicals. The start time of semiochemical usage was either on 19 February (50th day) or on 29 April (120th day) in the simulations for case studies ( Figure 6A). A variable duration of semiochemical release was generated following a normal distribution D~N(60, 20), where D represents duration in days ( Figure 6A, right). Since the dosage of semiochemicals was fixed (D·SCs = 720), the inner areas of each square were all the same (green, yellow to red, Figure 6A, right). Variations in the durations of release resulted in changes in the efficiency of attracting wasps or parameter λ ( Figure 6A, left). Stochastic simulations were then performed, and the relation between the tea yields and duration of semiochemical release was determined. Significant correlation was identified irrespective of the starting time ( Figure 6B,C, top). Increasing the duration of releases raised the tea yields in a year ( Figure 6B,C, bottom panels, p < 0.0001). Qualitatively similar results were obtained if semiochemicals were used as repellents ( Figure S4). Insects 2022, 13, x FOR PEER REVIEW 11 of 15  We further investigated whether there existed an optimum time for the application of semiochemicals. The starting time was varied from the 30th day to the 130th day. The averaged tea yields after using semiochemicals were then determined (Supplementary Materials: Average daily tea yields after semiochemical application). We found that the averaged tea yields first increased and then decreased when the starting time was changed from the 30th day to the 130th day ( Figure 7A-H). There was a local maximum around the 90th day ( Figure 7I). These results suggest that semiochemicals usage can be optimized to obtain higher average tea yields.

Discussion
In the current work, we constructed a model based on the tritrophic relationship among tea crops, TGLs, and parasitic wasps, and their interactions can lead to sustained oscillations. We further developed a realistic and stochastic model with temperaturedependent regulation. Generally, stochastic population models perform better than deterministic models in tritrophic ecosystems [43]. Our stochastic simulations showed that two leafhopper outbreaks occur in a year, which is compatible with field observations [44,45]. The first peak was relatively concise, whereas the second one was highly variable ( Figure 3C), which is consistent with previous reports [38,39]. Notably, diverse natural enemies such as entomopathogens, specialist and generalist predators coexist in tea ecosystems to exert biological control over tea pests [2]. Tritrophic relationship models can be regarded as interlocked predator-prey models. Even a minimal predator-prey model in tea ecosystems can exhibit complex dynamic patterns, for example, coexistence of multiple limit cycles and saddle node bifurcation [46]. Furthermore, multiple species of herbivores may also compete for a common resource (plant) [14]. Competition also occurs among natural enemies feeding on the same pests. Our model can be more precise if dynamic interactions among multiple predators or parasitoids and natural enemies are considered. Our model does not account for the immigration of leafhoppers, nor of their natural enemies, and these events may be influential on the dynamics of each population.
Previous models have usually ignored the effect of temperature dependence of prey or predator traits, especially in tritrophic interactions involving tea plants and TGLs [45][46][47]. In our current work, the temperature dependency of TGL developmental stages, lifespans, and tea growth rates were incorporated into the stochastic model. We found that effective accumulated temperature plays an important role in TGL outbreaks, with peaks identified within a specific range of average daily EATs ( Figure 5). Additionally, the distribution of frequency or possibility of pest outbreaks can be further decreased if the recording duration is increased and becomes relatively stabilized for n ≥ 80 ( Figure 5). Simulation results indicated that a long-term record of average daily EATs may be essential for predicting TGL outbreaks. If the long-term daily average EAT is above a certain threshold (e.g., approximately 18 °C in our model), the probability of pest outbreaks will be increased, and better pest management strategies should be implemented.
Sensitivity analysis showed that the systemic features are only sensitive to a few parameters ( Figure 1D). Specifically, r, K, and h represent the processes responsible for tea

Discussion
In the current work, we constructed a model based on the tritrophic relationship among tea crops, TGLs, and parasitic wasps, and their interactions can lead to sustained oscillations. We further developed a realistic and stochastic model with temperature-dependent regulation. Generally, stochastic population models perform better than deterministic models in tritrophic ecosystems [43]. Our stochastic simulations showed that two leafhopper outbreaks occur in a year, which is compatible with field observations [44,45]. The first peak was relatively concise, whereas the second one was highly variable ( Figure 3C), which is consistent with previous reports [38,39]. Notably, diverse natural enemies such as entomopathogens, specialist and generalist predators coexist in tea ecosystems to exert biological control over tea pests [2]. Tritrophic relationship models can be regarded as interlocked predator-prey models. Even a minimal predator-prey model in tea ecosystems can exhibit complex dynamic patterns, for example, coexistence of multiple limit cycles and saddle node bifurcation [46]. Furthermore, multiple species of herbivores may also compete for a common resource (plant) [14]. Competition also occurs among natural enemies feeding on the same pests. Our model can be more precise if dynamic interactions among multiple predators or parasitoids and natural enemies are considered. Our model does not account for the immigration of leafhoppers, nor of their natural enemies, and these events may be influential on the dynamics of each population.
Previous models have usually ignored the effect of temperature dependence of prey or predator traits, especially in tritrophic interactions involving tea plants and TGLs [45][46][47]. In our current work, the temperature dependency of TGL developmental stages, lifespans, and tea growth rates were incorporated into the stochastic model. We found that effective accumulated temperature plays an important role in TGL outbreaks, with peaks identified within a specific range of average daily EATs ( Figure 5). Additionally, the distribution of frequency or possibility of pest outbreaks can be further decreased if the recording duration is increased and becomes relatively stabilized for n ≥ 80 ( Figure 5). Simulation results indicated that a long-term record of average daily EATs may be essential for predicting TGL outbreaks. If the long-term daily average EAT is above a certain threshold (e.g., approximately 18 • C in our model), the probability of pest outbreaks will be increased, and better pest management strategies should be implemented.
Sensitivity analysis showed that the systemic features are only sensitive to a few parameters ( Figure 1D). Specifically, r, K, and h represent the processes responsible for tea plant growth and death or harvesting. Increasing the tea yields by promoting tea growth should not be a priority, as the equilibrium will shift toward higher pest populations via tritrophic relationships. Parameters µ (death rate of TGLs) and ω 1 (conversion rate from eggs to nymphs), were both sensitive to temperature [35]. However, the environmental temperature also influences various aspects of ontogenetic processes of leafhoppers and tea plant growth. As demonstrated above, it is not feasible to identify an optimum temperature to maximize tea growth while simultaneously attenuating leafhopper development with the highest efficiency. Furthermore, pesticide application can also decrease pest population (i.e., increasing µ) in the short term, but can subsequently result in a resurgence of pests, possibly with higher densities [48]. c 1 is a physiological attribute in TGLs and cannot be easily manipulated. The only possibilities are strategies targeting the insect feeding process (a 1 and b 1 ) and the attraction of natural enemies (λ). The application of synthetic semiochemicals seems to be a feasible way to increase tea yields in tritrophic systems. Repellants can directly decrease insect feeding, whereas attractants, for example, (E)-2-hexenal, benzaldehyde, and α-farnesene, can be used to lure natural enemies [9,10]. Mass-rearing and field release of parasitoid wasps also represent compelling strategies for integrated pest management by regulating λ [49]. All these strategies can indeed lead to a new balance with lower pest populations and provide better biocontrol performance.
With regard to the release of semiochemicals, our simulations suggested that sustained or slow release with low concentrations correlates with higher tea yields ( Figure 6). Slow-releasing semiochemicals outcompete fast-releasing regimes, irrespective of the initial application date. Our simulation might explain why volatile semiochemicals used as repellants or attractants should be coated with liquid paraffin to achieve better performance [11]. A possible mechanism is illustrated in Figure S5. When semiochemicals are released in short durations, disruption to pest population will be abated and pests will quickly revert to their original states (outbreak will still occur, Figure S5, red). When semiochemical application lasts longer, pest outbreak will be reduced, with increased tea crop biomass or tea yields ( Figure S5, blue). Collectively, using slow-releasing semiochemicals as repellants or attractants might be a compelling strategy in pest control.

Conclusions
In conclusion, our temperature-dependent stochastic model identified key features of the tritrophic relationship among tea plants, TGLs, and natural enemies in tea ecosystems. We further provided theoretical evidence that slow-releasing semiochemicals targeting insect feeding or the attraction of natural enemies may be highly efficacious. Using more temperature datasets and running the model for consecutive years may provide more insight into the dynamics of different populations in tea ecosystems under future global warming scenarios.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/insects13080686/s1. Figure S1: The curve-fitting to experimental measurements; Figure S2: The bifurcation diagram; Figure S3: The daily temperature at the year of 2020; Figure S4: Treatment with volatile semiochemicals to repel TGLs; Figure S5: Semiochemical application with different durations; Table S1: The parameter values and initial conditions; Table S2: Fitted parameters.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.