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 andPityohyphantes phrygianus) and twomirids (Nesidiocoris tenuis andMacrolophusmelanotoma). Biologicallymeaningful 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.


Introduction
A key challenge in ecology is to determine how changes in biodiversity affect ecosystem functioning and ecosystem service provision [1][2][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][5][6][7][8][9][10]. In fact, different aspects of enemy biodiversity including species identity [11][12][13][14], complementary prey use [6,[15][16][17], behavioural and ecological interactions between predators [18][19][20][21][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][25][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][31][32][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.

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 [36]. 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][38][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 [36].  [25] 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 [25]. 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 [29] or a response surface design [17] 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.

Experimental
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 [36], 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 [44]. 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 3 inspection. Insects were stored in alcohol and subsequently counted by means of a stereomicroscope. The number of whitefly adults ( Ad (O) ) 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 ( Ny-EU ). 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 ( Ny(O) ) 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.

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 [30].
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 ( Ny ) and adults ( Ad ) 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, (Ny) , and adults, (Ad) , 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][46][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 ( Ny-EU ) that survive predation on the nymph stage due to different predators: Ny-EU ∏ (1− (Ny) ) , 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 (∏ (1 − (Ad) ) ). 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 ( Ny-EU ) 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 27.7 ± 0.78 ∘ C; Block B 22.7 ± 1.07 ∘ C; = 6), 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 Ny-EU . 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.

ISRN Ecology
Finally, two parameters representing unknown mortality of the whitefly adults ( Ad ) and nymph subpopulations ( Ny ) 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.

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) [30].
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 ( Ny(O) ), whitefly adult emergence (Emerged), and whitefly adult survival ( Ad(O) , 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.

Single Series
Analysis. The first part of the processwhitefly nymph survival-was studied by analysing the series of counts Ny(O) , that is, the number of whitefly nymphs survived in leaf samples. These observations were predicted as follows: where Ny(E) 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 (+ NB and = 27, see below), was simplified according to the details presented in Table 1. 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 Emerged (E) 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 blockspecific emergence rates, the general model was simplified to specify no variation of emergence rates between blocks (⋅). The third part of the process-whitefly adult survivalwas studied by analysing the series of counts Ad(O) , 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.

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 [31].
The three series of counts ( Ny(O) + Emerged + Ad(O) ; = 66 * 3 = 198) 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.

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 [26]. 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 [48], and intraguild predation (IGP). Different means could arise from different predatorspecific detectability, different background mortality of predator species (e.g., reduction of predator numbers from 20 ISRN Ecology 5 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 (M IGP : IGP and different background mortality present), (b) average number of mirids recovered for each species (M BM : different background mortality present), and (c) overall average number of mirids recovered (M NULL : 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.

Error Structure and Model Selection.
The negative binomial distribution, suitable for analysis of counts [33,49,50] and invertebrate counts in particular [51], was used to specify the error structure in the data [31,33,52]. With the negative binomial distribution, the variance = mean + NB mean 2 . This distribution reduces to Poisson when NB = 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 [52] using the Newton search implemented in Microsoft Excel solver, which is a user-friendly and a reliable tool for maximum likelihood estimation [50]. 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 NB 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, QAIC , 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 / < 40) and when data are overdispersed [30,50]; where is the effective sample size, and log(L( | , )) 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 [52]. With this criterion, consistent departures from the value of 1 indicate inadequate fit. With real data, the estimated overdispersion parameter should generally be 1 ≥ ≤ 4 [30]. Overdispersion ( ), as well as the negative binomial NB , is estimated as a parameter and is therefore always considered in the , the number of estimated parameters of model . Moreover, the QAIC values were calculated by using a weighted estimate of overdispersion. This approach is fundamentally an application of model averaging [30]. 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 [30], = ∑̂.
Finally, support for each model was assessed by using ΔQAIC = QAIC ( ) − Lowest QAIC and the Akaike weight, , and a measure of the explanatory power ( 2 ) of the best models selected was estimated using the dispersion parameter NB of the following negative binomial [53]: where max is estimated in the negative binomial model with only a constant mortality term, ( ). Table 5  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 6 ISRN Ecology 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. Table 1 shows model selection results and parameter estimates coming from the analysis of data relative to whitefly nymph survival ( Ny(O) ; models S1-S17). The best model selected, S14 ( 2 = 0.9755):

What Do Models Capture?
(Ny) ( , , , , ) ; = = 0; specifies that there is no predation on nymphs from O. lineatus and P. phrygianus ( = = 0); there is some unknown mortality in all replicates ( ); per capita predation rates of N. tenuis and M. melanotoma vary with predator identity ( ); block effects are present and do not change with predator identity ( ); whitefly nymph starting density affects per capita predation rates differently with predator identity ( ); per capita predation rates differ between single-and multiplepredator treatments ( ). Specifically, when N. tenuis is in combination with any other predator species, per capita predation rates are higher than the single predator N. tenuis treatment ( < + ; = 0.1171, + = 0.1722). By contrast, M. melanotoma does not seem to change per capita predation rate when in combination with any other predator ( = + = 0.3338). However, model S13 ( = 0.4473; ΔQAIC = 0.0244) provides evidence for a positive interaction between M. melanotoma and P. phrygianus, as per capita predation rate of M. melanotoma increases when P. phrygianus is present ( = 0.2974; + = 0.4972). Figure 1 compares the number of whitefly nymphs counted in each leaf sample with expectations and standard deviations (SDs) provided by model S13 embedded in model J9 (see below, Table 2). The figure shows the range of variation present in the data and how the model is able to predict observed counts. As also suggested by the very high value of the estimated 2 , the model provides a good description of the variation present in the data.
The other models presented in Table 1 have substantially no support, either because they are overfitted or underfitted; that is, they contain too many or too few parameters [30]. However, models need to be presented in order to show supporting evidence for the effects retained in the best models. For example, simplification from model S14 to model S15 (ΔQAIC = 7.6905) provides strong support for the increase of N. tenuis predation rates in all multiple-predator treatments. To evaluate strength of evidence relative to prey density effects on predation rates ( ), see the simplification from S4 to S5. To evaluate strength of evidence relative to block effects on predation rates ( ), see the simplification from S14 to S16 (ΔQAIC = 25.7016). Also, the very high ΔQAIC of model S17 (57.8675), which assumes that there is only an unknown constant mortality for all replicates, provides strong overall support for the parameters retained in the best models.  Table 2 shows model selection results and parameter estimates coming from the joint analysis, with particular reference to whitefly adult emergence and predation on whitefly adults. According to the ΔQAIC , the information relative to predation on whitefly adults appears rather uncertain, and four models have similar weight (J8-J11). The best model selected, J9 ( 2 = 0.1037 and 2 = 0.6243, based on the NB of the Ad(O) and Emerged data series, resp.): (Ny) (S13) ( ) (Ad) ( ) ; specifies that predation on nymphs is as in model S13; emergence rate of whitefly adults changes with blocks, ( ), being less likely in Block B, where average day temperature was lower ( = 0.4443; = 0.2266); predation on whitefly adults varies between single-and multiple-predator treatments and does not appear to vary with predator identity, blocks, and starting whitefly nymph density, (Ad) ( ); there  The general model S1, (Ny) ( , , , , ) 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 → 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, NB , and overdispersion, ). Model selection diagnostics are described in the method section. Parameter estimates are: unknown mortality ( Ny ) 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.  (Ny) (S13) Ψ( ) (Ad) ( , , , , ) all, specifies that: model structure for whitefly nymph survival is as in S13, (Ny) (S13); emergence rates vary with blocks, Ψ( ); relatively to whitefly adult survival, (Ad) (⋅ ⋅ ⋅ ), 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 ( );

ISRN Ecology
parameter structure refers to all four predators, (all). The variation of model structure following each simplification step is: J1 →  In fact, the other three models sharing similar support state that (1) there is predation on whitefly adults from both spider species in single-predator treatments only ( = , J8); (2) spiders prey on whitefly adults in single-predator treatments and in combination, but not when they are with the mirids ( = = + = + , J10); there is no predation at all on whitefly adults, (Ad) ( ) (J11). Whereas models J8, J9, and J10 state similar things about spider predation on whitefly adults, it is important to note that ) and in combination with P. phrygianus ( + ) decrease with increasing whitefly nymph density. Per capita predation rates of N. tenuis alone ( ) and in combination with any other predator ( + ) tend to decrease less sharply with whitefly density. The two predators appear complementary across prey density. Parameter estimates taken from model J9. support obtained from model J11 and the very low 2 value suggest that predation on whitefly adults is a weak effect. Models J8, J9, and J10 also suggest that spiders further reduce their predation on whitefly adults to undetectable effects when predatory mirids are present. Support for the absence of substantial effects of mirids on whitefly adult survival is evident looking at simplification steps J3 to J4 and J5 (ΔQAIC that goes from 38.16 to 27.23 and 16.81, resp.).
On the other hand, it is important to note that the very high ΔQAIC of model J14, where emergence is assumed to be equal in Blocks A and B, provides strong evidence about variation of emergence rates of whitefly adults between temporal blocks. This is also confirmed by the 2 of the best model based on the NB relative to the Emerged data series, which is rather high, 0.6243, relative to model ( ) versus (⋅).

Functional Differences between the Predatory Miridae.
The well supported model parameters relative to mirid predation on whitefly nymphs provide some interesting information about possible functional differences on the two species investigated here. By using predation rates taken from model J9, Figure 3 shows how predation on whitefly nymphs by the predatory Miridae N. tenuis and M. melanotoma theoretically varies with whitefly nymph density and the presence of other predator species. Apparently, M. melanotoma is less effective than N. tenuis at high prey density and vice versa at low prey densities.

Predation between Predators.
As described in the methods section, all models ranked in Table 1 assume that number  of predatory Miridae changed during the experiment according to Table 4, thus, corresponding to the structure "M IGP : IGP and different background mortality present. " Model selection performed a posteriori to evaluate the presence of IGP between predators showed that, compared with model S14 (S14-M IGP ), model "S14-M BM : different background mortality present, no IGP" has a ΔQAIC = 5.5942, whereas model "S14-M NULL : no different mortality and no IGP" has a ΔQAIC = 5.5978. The rather high ΔQAIC of the two models supports the hypothesis that the numbers of mirid predators actually preying changed in each treatment as a consequence of different background mortality and IGP. However, considering that most of the ΔQAIC (5.5942) results from excluding IGP (S14-M BM ), the analysis suggests the presence of intraguild predation between some combinations of predators. Specifically, IGP between the predatory Miridae and from Linyphiidae towards N. tenuis is suggested.

Discussion
This study, performed in a relatively simplified experimental system, suggests several dimensions across which whitefly control is affected by natural enemies in tomatoes. We will briefly discuss below uncertainty and relevance of observed effects for pest control.

Different Predator Effectiveness.
In general, predatorspecific effectiveness related to identity and complementarity across prey life stages are well-known dimensions of predator biodiversity effects on prey [6,11,13,14]. Here, predation rates varied with predator identity and whitefly life stage, with the predatory Miridae and the spiders preferentially preying on whitefly nymphs and adults, respectively. Predatory Miridae have a strong and clear effect on whitefly nymphs. By contrast, spiders appear to have a weak effect on whitefly adults. Spiders are generalist predators known to feed infrequently and are often effective at low prey density [54,55], and their weak effect can be thus difficult to capture. However, the facts that predation from the Linyphiid spider P. phrygianus was directly observed during the experiment (dead whitefly adults were found trapped in spider webs present inside the experimental unit) and that both Linyphiidae and Oxyopidae are known to feed on small insects including Homoptera [40,55] support the evidence about predation from spiders on whitefly adults and suggest a role of spiders for whitefly control in tomatoes at low prey densities. [25,26,56] appeared as another dimension of variation of enemy biodiversity effects on whitefly control. The number of predatory Miridae recovered in multiple predator treatments was lower than that in single predator treatments, suggesting the presence of intraguild predation in some combinations. Also, interactions between predators could result in prey risk reduction and prey risk enhancement. Apart from IGP, evidence about prey risk reduction provided by model selection results (no predation of whitefly adults by spiders in the presence of predatory Miridae) is rather weak, and more evidence should be gained before discussion.

Contrasting Interactions between Predators. Ecological interactions between predators
On the opposite, prey risk enhancement is well supported by the constant increase of the estimated per capita predation rate of N. tenuis when combined with any of the other predators ( 2 and ΔQAIC strongly support this hypothesis; note also that this result holds with model S14-M NULL : no different mortality and no IGP, where the density of predators is constant). Variation of within-plant distribution of mirid predation on whitefly nymphs in multiple-predator treatments due to predator interference was shown in a similar system with direct observations [57]. Higher mobility of predators and exploitation of dispersed whitefly nymph patches could result in increased complementarity [26] and explain prey risk enhancement. However, as we cannot rest on direct observations, this mechanism should be seen as a hypothesis to be further evaluated with more focused auxiliary observations.

Functional Differences between Predatory Miridae.
In general, knowledge about functional response curves of predators is important in order to properly describe multiple predator effects [58,59]. Here, whitefly nymph density differently affected per capita predation rates of predatory Miridae ( Figure 3). Specifically, M. melanotoma appears to be less effective than N. tenuis at high prey density, and vice versa at low prey density. The high handling time and the poor capability to aggregate in high prey density patches of M. melanotoma (caliginosus) [60] and the enhanced plant feeding of N. tenuis with low whitefly density [38] may explain this pattern. The two mirid species thus appear complementary along a prey density axis and appear to have different roles in whitefly  Table 2, bottom), per capita predation rates on whitefly nymphs of the two predatory Miridae are lower in Block A. The effect found might be due to (a) increased phytophagy, which has been shown in N. tenuis to be related with high temperatures [38] and (b) prey reaching a refuge from mirid predation by transiting from the nymph to the adult stage faster. Note that the likely negative effect on per capita predation rates of the overall higher whitefly nymph density observed in Block A should have been already captured by the scaling factor . Given the utility of refuges for prey survival (e.g., [26]), environmental temperatures might contribute to whitefly outbreaks in warm conditions [63] by differentially affecting population growth of whiteflies (nymphs "refuging" in the adult stage) and natural enemy performance (more phytophagy and lower predation).

Enemy Biodiversity and Whitefly
Control. From a general enemy biodiversity and pest control perspective, this study suggests that the predator species evaluated here have complementary roles in whitefly control. The two species of predatory Miridae are both very effective on whitefly nymphs but showed different functional responses and effectiveness with different prey densities. Also, the mirids seem unable to prey on whitefly adults. The spiders, by weakly preying on adults, may contribute a small but important component of total whitefly suppression, with particular reference to the beginning of the season when whiteflies colonize crop fields and prey densities are low. Moreover, as multiple predator treatments did not show a reduction of whitefly control, spiders do not appear to interfere considerably with the predatory Miridae. A reduced relevance of intraguild predation for mirid population growth and whitefly control is also confirmed by the very high densities generally reached by predatory Miridae in untreated and enemy-rich tomatoes [36]. In contrast to other systems (e.g., [18,64]), we have no evidence to state here that realistic densities of intraguild predators disrupt biological control [21,65]. Therefore, even from a limited perspective of singleprey suppression and short term experiment, and by considering only pair-wise interactions, assembly of the four predator species appears here rather like an accumulation of component parts of biological control function, with predator complementarity arising across prey life stages and density and the emergence of weak interactions between predators that did not result in whitefly control disruption. These results suggest that increasing enemy biodiversity and whitefly control in tomatoes can be compatible and related goals [66]. At the same time, it must be acknowledged that interaction strength between natural enemies may greatly vary in different contexts (e.g., [24]). In fact, the net result of the interaction between natural enemies and prey in agroecosystems has always some level of uncertainty and context dependence. This is an intrinsic problem of enemy biodiversity research and pest control applications, which can be dealt with by properly designing and monitoring sustainable pest management [67].

Conclusions
By providing a structural description of the roles that different natural enemies have on whitefly suppression, this study informs us about biological control of whiteflies in tomatoes. As the pest control functioning of agroecosystems is affected by locally different biological communities and ecological conditions, site-or regional-specific information about agroecosystem composition, structure, and functioning is needed to foster sustainable agriculture [27,68]. In fact, pest management in crop systems goes on despite considerable lack of knowledge and uncertainty about agroecosystem functioning [67]. Decisions are made on a daily basis to guarantee food production, and available information about agroecosystem functioning should be used at best to drive hypotheses about sustainable pest management options. Within this applied perspective, the confrontation of realistic process models with available data can be useful to evaluate and develop hypotheses about the local structure of enemy biodiversity effects on prey, providing information about pest control functioning. To estimate biologically meaningful parameters, auxiliary information [33], such as stage transitions, predation on different stages, spatial and behavioural patterns of predation, variation of predator and prey densities and environmental conditions, can be measured and embedded in process models. Well supported parameters can then be used to make predictions and develop pest management scenarios to further evaluate hypotheses about local agroecosystem functioning by confronting management results with expectations [67].