Research Article | Open Access
Evaluating the Structure of Enemy Biodiversity Effects on Prey Informs Pest Management
Evaluating the structure of enemy biodiversity effects on prey in agroecosystems can provide insights into biological pest control functioning. With this aim, theoretical models that describe biological mechanisms underlying prey suppression can be developed and confronted with experimental data by means of model selection. Here, we confront multiplicative risk models to evaluate the structure of multiple predator effects on the whitefly Trialeurodes vaporariorum provided in tomatoes by two spiders (Oxyopes lineatus and Pityohyphantes phrygianus) and two mirids (Nesidiocoris tenuis and Macrolophus melanotoma). Biologically meaningful parameters retained in the best models showed that several predator traits differently affected pest control: species-specific per capita predation rates, prey use extent, different type of interactions between predators, and the response of predator species to prey density and environmental temperature. Even from a limited perspective of single-pest control and short term experiment, this study suggests that assembly of the four predator species results in predator complementarity across prey life stages and density, interactions of prey and predators with environmental conditions, and interactions between predators that do not result in whitefly control disruption. Such information about enemy biodiversity and whitefly control functioning can drive hypotheses about sustainable pest management options in local agroecosystems.
A key challenge in ecology is to determine how changes in biodiversity affect ecosystem functioning and ecosystem service provision [1–3]. With respect to the role of natural enemy biodiversity on pest control, biodiversity and ecosystem function studies typically explore the extent to which prey suppression is affected by variation of enemy biodiversity [4–10]. In fact, different aspects of enemy biodiversity including species identity [11–14], complementary prey use [6, 15–17], behavioural and ecological interactions between predators [18–22], and abundance [10, 23, 24] can determine emergent effects of enemy biodiversity on prey. Moreover, it is also becoming increasingly evident that emergent properties of enemy communities are context dependent and often arise from several simultaneous mechanisms that individually might produce linear or nonlinear, negative, null, and positive effects on prey suppression [5, 8, 9, 24–26]. It is therefore important to investigate multiple mechanisms underlying the link between enemy biodiversity and prey suppression to parse out the main processes determining pest control service provision in diverse agroecosystems [17, 26, 27].
To date, attempts to address the complexity underlying enemy biodiversity effects on prey have led to varied experimental approaches to partition individual mechanisms underlying prey suppression (e.g., [5, 8, 17, 28, 29]). Nevertheless, intrinsic features of living systems and experimental limitations make it difficult and often impossible to fully disentangle the simultaneous aspects of biological diversity underlying enemy biodiversity effects, resulting in multiple working hypotheses regarding experimental outcomes. To address this problem of multiple mechanisms, theoretical models that structurally attempt to describe the suite of possible processes involved can be developed, and relative support for hypotheses can be evaluated by confronting models with data [30–33]. A model selection approach based on realistic process models aims at developing relatively good approximations of the structure of biological processes [34, 35]. By focusing on the estimation of biologically meaningful parameters [31, 35], such approach can be a powerful tool to capture the information about the structure and magnitude of enemy biodiversity effects on prey and make predictions about system functioning to inform pest management.
With this aim, we extend and adapt multiplicative risk models [5, 25, 26] to provide a realistic structural description of predator biodiversity effects on the suppression of the whitefly Trialeurodes vaporariorum in tomatoes. The study focuses on four predator species that are common in tomato agroecosystems in south-eastern Sardinia (Italy): the actively hunting spider Oxyopes lineatus, the web-building spider Pityohyphantes phrygianus, and the predatory Miridae Nesidiocoris tenuis and Macrolophus melanotoma. Hypotheses about system functioning, that is, models, were informed by a priori understanding of predator biology and interactions with a focus on the main biological traits of predators thought to affect whitefly control. Model selection and maximum likelihood parameter estimates showed that several predator traits simultaneously affected whitefly control.
Overall, this paper attempts to show how process modeling can be useful to capture information contained in available data and inform pest management with a set of plausible hypotheses about pest control functioning in given local systems.
2.1. Biological System
Experiments were performed during Summer 2006 in a fresh-market tomato crop (variety Empire, Petoseed) situated in south-east Sardinia, Italy (NW corner — UTM: 32S 0524770; 4349400). In this area, the greenhouse whitefly Trialeurodes vaporariorum Westwood (Homoptera: Aleyrodidae) is often the most abundant pest . Two mirid predators, Nesidiocoris tenuis Reuter and Macrolophus melanotoma Costa (Heteroptera), feeding mainly on whitefly nymphs, are known to play a key role in whitefly control in Mediterranean tomatoes [37–39]. Also, the two spiders Oxyopes lineatus Latreille (Araneae: Oxyopidae) and Pityohyphantes phrygianus Koch (Araneae: Linyphiidae), by preferentially preying on mobile animals [40, 41], can prey on whitefly adults and can also act as intraguild predators by preying on predatory Miridae [42, 43]. Consequently, spiders potentially can both contribute to pest control with direct predation on whitefly adults and interfere with it by interacting antagonistically with the predatory Miridae. As these generalist predators are very common in south-eastern Sardinia open-field tomatoes, they likely represent an important part of local biodiversity affecting whitefly control .
2.2. Experimental Design
The effect of single-predator (N. tenuis, M. melanotoma, O. lineatus, and P. phrygianus; thereafter, Nt, Mm, Ox, and Li, resp.) and two-predator (Nt + Mm, Nt + Ox, Nt + Li, Mm + Ox, Mm + Li, and Ox + Li) treatments on whitefly survival was studied. Each treatment and a control (no predators) were replicated 3 times per 2 temporal blocks (A and B) for a total of 66 experimental units. An additive design  was used, where intraspecific density is kept constant in multiple predator treatments, whereas interspecific density increases. This design was selected to avoid confounding effects of variation of intraspecific densities, on the basis that spiders and mirids are functionally different and should not be considered substitutable . Nevertheless, it must be underlined that in any biodiversity experiment, changes in diversity necessarily produce changes in species densities (intra- or interspecific). Therefore, in order to fully address confounding effects of predator density variation, both the additive and the substitutive designs  or a response surface design  should be used. As in many experiments, this was not possible here for economic constraints. Nevertheless, as ecological data on agroecosystem functioning are rare and costly, it is always worth attempting to capture the information contained.
The experimental field was settled and managed by a local farmer (Claudio Orrù) in his own farm under the supervision of the authors. Two thousand tomato seedlings were transplanted at the beginning of June. Approximately one month later, 66 experimental units (tomato plants) were randomly selected and caged in mesh net-bags (0.5 mm mesh, 100 cm tall, and 40 cm in diameter) secured to the ground. Before cage closure, plants were treated with pyrethrum and visually inspected to eliminate arthropods. Five days after the closure of mesh net-bags, about 80 whitefly adults were introduced in each cage to start the whitefly population with egg deposition. Whitefly adults were collected 1 day before their release from an organic vegetable garden, managed by Paolo Casula and situated a few kilometers away, and stored at 5°C. At the beginning of whitefly adult emergence, which started approximately 20 days after whitefly introduction and egg deposition (July 25th and August 2nd for blocks A and B, resp.), predators were released and kept in experimental units for 6 days. In this way, predation on both the whitefly nymph and the adult stage could be studied.
To keep starting densities of single-and-multiple predator treatments within ranges commonly found in this ecological system , 20 individuals each of the predatory Miridae Nt and Mm (3rd instar nymphs), 5 (Li), and 4 (Ox) immature spiders of similar sizes were introduced in the experimental units. These densities can be observed in the study system during May-June, at the beginning of the season. As whiteflies colonize the crop early in the season, this is a key moment in natural pest control. Field observations show that later on in the season, mirids become the dominant predator species, reaching very high proportions with respect to any other natural enemy (more than 30 mirids per spider, mostly Nesidiocoris tenuis). Individuals of Nt were collected from a tomato field established within the same farm on early May, whereas Mm, Ox, and Li were collected in nearby seminatural habitats from Dittrichia viscosa (Asteraceae). This plant is very common in abandoned agricultural fields and is also known to host M. melanotoma . All predators were collected 1–3 days before release and stored at 5°C.
At the end of the experiment, cages were opened, and insects recovered from plants by means of a suction sampler (Vortis, Burkard Manufacturing Ltd.) and careful visual inspection. Insects were stored in alcohol and subsequently counted by means of a stereomicroscope. The number of whitefly adults () and predators (, , , and ) detected at the end of the experiment in each experimental unit was recorded (Table 5). All leaves were also visually inspected to count the total number of whitefly nymphs (dead, alive, and pupal cases from which whitefly adults emerged) present in the experimental unit at the end of the experiment (). This count was used as a best estimate of the total number of alive whitefly nymphs present in each experimental unit at the beginning of the experiment.
Additionally, samples of up to 12 infested leaves per experimental unit were observed under a stereomicroscope to count the number of nymphs dead and alive and pupal cases from which adult whiteflies emerged (Dead, Alive, and Emerged, Table 5). The number of survived nymphs observed in leaf samples () is then obtained by adding Alive with Emerged. This additional information, derived from leaf samples, was used to estimate parameters relative to predation on whitefly nymphs and adult emergence.
2.3. Modeling Enemy Biodiversity Effects on Prey
The starting point of model selection approaches is to develop a set of a priori hypotheses about the biological system expressed as mathematical or statistical models, so that likelihood-based approaches can be used for confrontation of models with data [30, 31, 33]. By drastically reducing the possibility to find spurious correlations, a priori modeling is fundamental to avoid data dredging .
When dealing with predation, attempts to develop realistic process models can start from multiplicative risk framework [25, 26], which needs to be adapted considering the specific study system and experimental setting. As described before, the experiment started when whitefly adult emergence began, and the number of whitefly nymphs () and adults () present in each experimental unit could be approximated by the following: During the experiment, the relative numbers of the two whitefly subpopulations change as a consequence of emergence rates, , and predation. Predation results in a reduction of the two subpopulations proportional to predator density, , and predator-specific per capita predation rates on whitefly nymphs, , and adults, , where is the predator species.
If the number of whiteflies surviving predation from a given predator species is approximated by the following general multiplicative risk model [25, 26, 45–47]: the number of whitefly adults expected at the end of the experiment can be written as follows: This model means that the number of adult whiteflies expected at the end of the experiment is given by the number of whitefly nymphs present at the beginning of the experiment () that survive predation on the nymph stage due to different predators: , where is predator species , , , or . Surviving nymphs become adults with rate , and to be counted at the end of the experiment, they have to survive predation on the adult stage (). By specifying in the general model above the relative density of each predator present in the experimental unit, specific cases of the model valid for each experimental unit can be obtained. With this specification, the general model has 10 estimated parameters (8 predation rates and 2 block-specific emergence rates).
As often happens in field studies, not everything could be controlled; that is, starting whitefly nymph densities (see Table 5) and environmental conditions of the two temporal blocks varied. On the other hand, measured variation might contain interesting information. The variation of the number of whitefly nymphs present in each experimental unit () was thus considered as a covariate in the analysis, allowing us to study the effect of starting prey density on predator-specific per capita predation rates. Also, the presence of strong cold wind greatly lowered the average environmental temperature during the experiment in block B (average temperature ± SD: Block A °C; Block B °C; ), possibly affecting whitefly adult emergence and predation. Therefore, the general model was further developed to consider likely effects of blocks (temperature) and starting whitefly nymph density on per capita predation rates as follows: where is an additive factor that can vary with predator identity and applies only to block A, whereas is a scaling factor that can vary with predator identity and relates per capita predation rate with whitefly nymph density . With this specification, the general model now has 2 more estimated parameters for each of the 8 predation rates: 8 + 2 + 16 = 26.
In order to model possible variations of per capita predation rate due to interactions between predators, per capita predation rates of each predator species when released in single- or in two-predator treatments were estimated separately as follows: where refers to the presence of another predator . With this specification, the general model now has 3 more estimated parameters for each predator species and prey stage: 8 + 2 + 16 + 24 = 50.
Finally, two parameters representing unknown mortality of the whitefly adults () and nymph subpopulations () can also be specified (52 parameters in the general model). Unknown mortality, which does not depend on the presence of predators, is assumed to be present in all experimental units. The 52 parameters contained in the full a priori model now need to be supported by the data to be kept for inference, as described below.
2.4. Model Simplification and Selection
In the previous section, we developed a general model embedding the main factors thought to affect variation in the data. During confrontation with data, the general model is simplified to search for parsimonious models actually supported by the data, achieved here by means of model selection based on Akaike’s Information Criterion (AIC) .
Starting from the general model described above, simplified models were developed by constraining or eliminating a single factor at a time. Factors were excluded if model simplification resulted in lower AIC values. If AIC values were higher, similar, or with lower than 1, factors were kept and evaluated again later on. For simplicity, the analysis was done separately on three series of counts relative to whitefly nymph survival (), whitefly adult emergence (Emerged), and whitefly adult survival (, see Table 5). The three data series described different parts of the overall whitefly survival process. The structure of the joint model was subsequently evaluated by reanalysing the three data series jointly.
2.4.1. Single Series Analysis
The first part of the process—whitefly nymph survival—was studied by analysing the series of counts , that is, the number of whitefly nymphs survived in leaf samples. These observations were predicted as follows: where is the expected number of whitefly nymphs present in leaf samples that survived at the end of the experiment and is the total number of nymphs counted in the sample (see Table 5). Predation rates were estimated as in models 4 and 5 developed in the previous section. The general model relative to this part of the process, which has 25 parameters and = 27, see below), was simplified according to the details presented in Table 1.
|The general model S1, all, specifies that: there is a constant unknown mortality of whitefly nymphs in all replicates (); predation rates on whitefly nymphs varies with predator identity (), blocks ( indicating that varies with predator identity) and whitefly nymph density ( indicating that varies with predator identity); species-specific predation rates vary also between single- and multiple-predator treatment (); parameter structure refers to all four predators, (all). The variation of model structure following each simplification step is: S1 → S2: no predation from O. lineatus is assumed (all → ); S2 → S3: no predation from P. phrygianus is further assumed (); S3 → S4: block effect constant between predator species ( → B); S4 → S5: effect of starting whitefly density on predation rates constant between predators (); S4 → S6: unknown mortality absent ( → ); S4 → S7: per capita predation rates of M. melanotoma are identical in Mm and Mm + Nt treatments (), and so on (S7 → S8 → S9; S8 → S10; S8 → S11; S8 → S12; S8 → S13 → S14 → S15; S14 → S16). : number of model parameters (including the dispersion parameter of the negative binomial distribution, , and overdispersion, ). Model selection diagnostics are described in the method section. Parameter estimates are: unknown mortality () and per capita predation rates of predators in single-(e.g., for M. melanotoma) and multiple-predator treatments (e.g., for M. melanotoma when in combination with N. tenuis). The sign — indicates that the parameter has been set to zero. Best models in bold. |
The second part of the process—whitefly adult emergence—was studied by analysing the series of counts Emerged, that is, the number of whitefly adults emerged in leaf samples. These observations were predicted as follows: where is the expected number of whitefly adults emerged in the sample, s is the total number of nymphs counted in the sample, and indicates block-specific emergence rates (2 parameters). To evaluate support for block-specific emergence rates, the general model was simplified to specify no variation of emergence rates between blocks .
The third part of the process—whitefly adult survival—was studied by analysing the series of counts , that is, the number of whitefly adults survived and recovered in the whole experimental unit. These observations were predicted using general model 3 specified in the previous section. The general model, which has identical structure and number of parameters as model 6 above relative to whitefly nymph survival, was subsequently simplified according to the details presented in Table 2. To run this model, parameter estimates relative to predation on whitefly nymphs and adult emergence were taken from previous single-series analysis.
|The general model J1, (S13) all, specifies that: model structure for whitefly nymph survival is as in S13, (S13); emergence rates vary with blocks, ; relatively to whitefly adult survival, , there is a constant unknown mortality in all replicates ; predation varies with predator identity , blocks ( indicating that varies with predator identity) and whitefly density ( indicating that varies with predator identity); species-specific predation rates vary also between single- and multiple-predator treatment ; parameter structure refers to all four predators, (all). The variation of model structure following each simplification step is: J1 → J2: unknown mortality absent ( → ); J2 → J3: block and starting whitefly density effects absent ( → ); J3 → J4: no predation from N. tenuis is assumed (all → ); J4 → J5: no predation from M. melanotoma is further assumed (→), and so on (J5 → J6; J5 → J7; J5 → J8; J5 → J9; J5 → J10; J9 → J12; J9 → J13). : number of model parameters (including the dispersion parameter of the negative binomial distribution, , and overdispersion, ). Model selection diagnostics are described in the method section. Parameter estimates presented are: emergence rates for block A and B ( and ) and per capita predation rates of spiders in single- (e.g., for O. lineatus) and multiple-predator treatments (e.g., for O. lineatus when in combination with P. phrygianus). At the bottom of the table is also presented the full set of parameter estimates from model J9 relative to predation on whitefly nymphs. The sign — indicates that the parameter has been set to zero. Best models in bold.|
2.4.2. Joint Analysis
Events resulting in whitefly nymph mortality, adult emergence, and adult mortality can be seen as a sequence of independent processes acting simultaneously. In these cases, likelihood theory can incorporate multiple sources of observations; that is, the likelihood of the joint model is the product of the likelihood of each series of data .
The three series of counts ( + + ; ) were thus confronted with expectations by using the joint likelihood of model 3, which embeds models 4, 5, 6, and 7. Estimates of parameters from the joint models were obtained by maximizing the sum log-likelihood of the three series of counts. The assumption of independence of predation on whitefly nymphs, emergence from pupae, and predation on adults is an approximation that allows an easier statistical analysis. The analysis performed on single and multiple sources of data showed similar model structure and parameter estimates, indicating that the assumption of independence did not affect results and was therefore not adjusted.
2.4.3. Evaluating Predation between Predators
There is another source of variation possibly affecting the process and the experimental outcomes: the number of predators actually preying on whiteflies, which could be reduced by intraguild predation. To address this possible variation, the exponents of general model 3 can be interpreted as known variables measured by the average number of predators recovered at the end of the experiment for each species and treatment (Table 4). This auxiliary information can be embedded in models to specify hypotheses about the real number of predators active during the experiment . In fact, the average number of predators detected seems to change with predator identity and treatment. If averages presented in Table 4 contain useful information, model likelihood should improve when specified in the models as exponents .
The low number of predators recovered at the end of the experiment could be due to predator background mortality, detectability lower than 1 , and intraguild predation (IGP). Different means could arise from different predator-specific detectability, different background mortality of predator species (e.g., reduction of predator numbers from 20 released to on average 8.8 and 6.2 recovered for and , resp.), and the presence of IGP in some combinations (e.g., reduction of predator numbers in multiple predator treatments from 8.8 and 6.2 to 5.5 and 3.5 for and , resp.).
The average numbers of predators recovered at the end of the experiment for each species and treatment were thus embedded in all models presented in Tables 1 and 2. The average numbers presented in Table 4 were approximated to the closer integer for predatory Miridae, whereas spider density was set to 1. The models therefore assume a variation of mirid density due to IGP and differential apparent background mortality between predators (detectability cannot be estimated and is therefore ignored, likely resulting in overestimation of per capita predation rates). Support relative to the presence of differential background mortality and IGP was evaluated by confronting the best model selected from the a priori set with models that specify as exponents: (a) average number of mirids recovered in single versus multiple treatments for each species (MIGP: IGP and different background mortality present), (b) average number of mirids recovered for each species (MBM: different background mortality present), and (c) overall average number of mirids recovered (MNULL: no different mortality and no IGP). As results relative to whitefly adult survival were very uncertain, this analysis was performed on the nymph survival data series only.
2.4.4. Error Structure and Model Selection
The negative binomial distribution, suitable for analysis of counts [33, 49, 50] and invertebrate counts in particular , was used to specify the error structure in the data [31, 33, 52]. With the negative binomial distribution, the variance = mean + mean2. This distribution reduces to Poisson when = 0. This estimate of variance was used to compute standard deviation of expected counts.
Maximum likelihood estimation of model parameters was achieved by numerical search  using the Newton search implemented in Microsoft Excel solver, which is a user-friendly and a reliable tool for maximum likelihood estimation . In practice, the likelihood of the model is given by the sum of the likelihood of each observation relative to model expectations, with the distribution of data points being assumed to follow the negative binomial with estimates of specific for each data series.
To address for possible overdispersion, very frequent in real insect count data, model selection was achieved through the calculation of quasi-likelihood adjustments for overdispersion () of the corrected Akaike’s Information Criterion, , which is a modified version of the AIC used when the number of the estimable parameters in the model () is large relative to sample size (i.e., when ) and when data are overdispersed [30, 50]; where is the effective sample size, and is the log-likelihood evaluated at the maximum likelihood estimates, , given data and model . Overdispersion is estimated as follows: where is sample size and is the number of estimated parameters of model . With this criterion, consistent departures from the value of 1 indicate inadequate fit. With real data, the estimated overdispersion parameter should generally be . Overdispersion (), as well as the negative binomial , is estimated as a parameter and is therefore always considered in the , the number of estimated parameters of model . Moreover, the values were calculated by using a weighted estimate of overdispersion. This approach is fundamentally an application of model averaging . As overdispersion is an estimated parameter clearly affected by model structure, it seems sounder to estimate the parameter by taking into account model weight. The weighted average of overdispersion was thus obtained by using Akaike weights , .
Finally, support for each model was assessed by using and the Akaike weight, , and a measure of the explanatory power () of the best models selected was estimated using the dispersion parameter of the following negative binomial : where is estimated in the negative binomial model with only a constant mortality term, .
3.1. What Do Raw Data Suggest?
Table 5 shows how counts of the 3 data series analysed varied among experimental units and between blocks: = (SD) and for blocks A and B, respectively; Emerged = and for A and B; = and for A and B. The range of variation observed is a consequence of both the great variation of the number of whitefly nymphs present in each experimental unit at the beginning of the experiment ( = and ; ranges 213–1776 and 30–1230 for Blocks A and B, resp.), likely due to heterogeneity in whitefly population growth after adult introduction in cages (80 per cage), and block effects. In fact, mean numbers observed and some simple proportions suggest that emergence rates (Emer./s (tot) = 0.38 and 0.23 and / = 0.40 and 0.14 for Blocks A and B, resp.) and nymph mortality (Dead/s (tot) = 0.16 and 0.28) differ between blocks, suggesting slower adult emergence and higher predation in Block B.
Mean mortality rates of whitefly nymphs estimated for each treatment are shown in Table 3 (Dead/s (tot), Table 5 footer). The table suggests that the two predatory Miridae N. tenuis (Nt) and M. melanotoma (Mm) have different effects, whereas the effect of spiders does not seem to differ with the control. Interestingly, whitefly nymph mortality observed in predator combinations seems determined by the most effective predator and is not increased or drastically reduced by the presence of other predators. However, raw data presentation cannot capture all the information contained, and the model selection results will show below a more complex picture.
|Dead: number of dead whitefly nymphs counted in leaf samples; : sample size (Dead + Alive + Emerged); Nt, Mm, Ox, Li: predator treatments (N. tenuis, M. melanotoma, O. lineatus, and P. phrygianus resp.), C: control. |
|Number of predators recovered for each experimental unit are given in Table 5 (, , , ), where Nt, Mm, Ox, Li are predator treatments (N. tenuis, M. melanotoma, O. lineatus, and P. phrygianus resp.); Mean: average number of predators recovered for each species.|
|: number of whitefly nymphs counted in each experimental unit; : number of whitefly adults counted in each experimental unit; , , ,|