# Evaluating the Structure of Enemy Biodiversity Effects on Prey Informs Pest Management

Authors: Paolo Casula; Mauro Nannini
Journal: ISRN Ecology
 (2013)

---

## Abstract

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 whiteflyTrialeurodes 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.


---

## Body

## 1. Introduction

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. Methods

### 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 whiteflyTrialeurodes 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–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].


### 2.2. Experimental Design

The effect of single-predator (N. tenuis, M. melanotoma, O. lineatus,and P. phrygianus;thereafter, Nt, Mm, Ox, andLi, resp.) and two-predator (Nt + Mm, Nt + Ox, Nt + Li, Mm + Ox, Mm + Li,andOx + 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 [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.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 inspection. Insects were stored in alcohol and subsequently counted by means of a stereomicroscope. The number of whitefly adults (NAd(O)) and predators (NNt, NMm, NLi, and NOx) 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 (NNy-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 (NNy(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.


### 2.3. Modeling Enemy Biodiversity Effects on Prey

The starting point of model selection approaches is to develop a set ofa 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 (NNy) and adults (NAd) present in each experimental unit could be approximated by the following:
(1)NNy≈NNy-EU,NAd≈0.During the experiment, the relative numbers of the two whitefly subpopulations change as a consequence of emergence rates, ψ(B), and predation. Predation results in a reduction of the two subpopulations proportional to predator density, Ni, and predator-specific per capita predation rates on whitefly nymphs, pi(Ny)′, and adults, pi(Ad)′, where i is the predator species.If the number of whiteflies surviving predation from a given predator speciesi is approximated by the following general multiplicative risk model [25, 26, 45–47]:
(2)NWhitefly(t+1)=NWhitefly(t)(1-pi′)Ni,
the number of whitefly adults expected at the end of the experiment can be written as follows:
(3)NAd(E)=NNy-EU(1-pNt(Ny)′)NNt(1-pMm(Ny)′)NMm×(1-pOx(Ny)′)NOx(1-pLi(Ny)′)NLiψ(B)×(1-pNt(Ad)′)NNt(1-pMm(Ad)′)Nm×(1-pOx(Ad)′)NOx(1-pLi(Ad)′)NLi.
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 (NNy-EU) that survive predation on the nymph stage due to different predators: NNy-EU∏in‍(1-pi(Ny)′)Ni, where i is predator species Nt, Mm, Ox, or Li. Surviving nymphs become adults with rate ψ(B), and to be counted at the end of the experiment, they have to survive predation on the adult stage (∏in‍(1-pi(Ad)′)Ni). 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 Table5) 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 (NNy-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; n=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:
(4)pi(block,unit)′=pi′+βie(NNy-EU*ai),
where βi is an additive factor that can vary with predator identity and applies only to block A, whereas ai is a scaling factor that can vary with predator identity and relates per capita predation rate with whitefly nymph density NNy-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:(5)pi(Ny)′≠pi+Pj(Ny)′,pi(Ad)′≠pi+Pj(Ad)′,
where Pj refers to the presence of another predator j. 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 (pCAd) and nymph subpopulations (pCNy) 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) [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 (NNy(O)), whitefly adult emergence (Emerged), and whitefly adult survival (NAd(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.


#### 2.4.1. Single Series Analysis

The first part of the process—whitefly nymph survival—was studied by analysing the series of countsNNy(O), that is, the number of whitefly nymphs survived in leaf samples. These observations were predicted as follows:
(6)NNy(E)=s(1-pNt(Ny)′)NNt(1-pMm(Ny)′)NMm×(1-pOx(Ny)′)NOx(1-pLi(Ny)′)NLi(1-pCNy),
where NNy(E) is the expected number of whitefly nymphs present in leaf samples that survived at the end of the experiment and s 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 (+kNB and c = 27, see below), was simplified according to the details presented in Table 1.Table 1
Best model structure and main parameter estimates relative to enemy biodiversity effects on whitefly nymph survival.
Code
Model structure
Model selection diagnostics
Maximum likelihood estimates of model parameters
ΔQAICc
w
i
K
i
c
i
R
2
i
k
NB
p
C
Ny
p
M
m
′
p
M
m
+
N
t
′
p
M
m
+
O
x
′
p
M
m
+
L
i
′
p
N
t
′
p
N
t
+
O
x
′
p
N
t
+
L
i
′
p
N
t
+
M
m
′
S1
p
(
Ny
)
(
C
,
i
,
B
i
,
a
i
,
M
);  all
57.6253
0.0000
27
2.3675
1.0000
0.0000
0.0729
0.3004
0.0000
0.0000
0.4914
0.1611
0.0256
0.1192
0.2908
S2
p
(
Ny
)
(
C
,
i
,
B
i
,
a
i
,
M
)
;
p
O
x
′
=
0
34.4510
0.0000
21
2.2286
0.9944
0.0005
0.0869
0.2860
0.4184
0.2474
0.4899
0.1097
0.1541
0.1787
0.1552
S3
p
(
Ny
)
(
C
,
i
,
B
i
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
14.2799
0.0004
15
1.8687
0.9810
0.0015
0.0894
0.3128
0.3905
0.2870
0.5029
0.1163
0.1633
0.1802
0.1700
S4
p
(
Ny
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
11.5475
0.0014
14
1.8225
0.9787
0.0017
0.0914
0.3001
0.3798
0.2615
0.4784
0.1148
0.1614
0.1785
0.1678
S5
p
(
Ny
)
(
C
,
i
,
B
,
a
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
26.3280
0.0000
13
1.6959
0.9072
0.0075
0.0906
0.1021
0.1645
0.0850
0.1304
0.1336
0.2037
0.2306
0.1646
S6
p
(
Ny
)
(
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
24.1824
0.0000
13
1.4922
0.8872
0.0092
—
0.2716
0.2591
0.2200
0.3759
0.1208
0.1598
0.1878
0.1662
S7
p
(
Ny
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
N
t
′
8.3366
0.0070
13
1.7840
0.9783
0.0018
0.0914
0.3019
=pMm′
0.2616
0.4788
0.1159
0.1633
0.1804
0.1739
S8
p
(
Ny
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
N
t
′
=
p
M
m
+
O
x
′
5.3157
0.0317
12
1.7565
0.9785
0.0017
0.0917
0.2976
=pMm′
=pMm′
0.5006
0.1162
0.1638
0.1809
0.1757
S9
p
(
Ny
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
P
j
′
5.1225
0.0350
11
1.7795
0.9759
0.0020
0.0917
0.3323
=pMm′
=pMm′
=pMm′
0.1152
0.1619
0.1792
0.1704
S10
p
(
Ny
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
N
t
′
=
p
M
m
+
O
x
′
;
p
N
t
′
=
p
N
t
+
O
x
′
8.1180
0.0078
11
1.6109
0.9515
0.0039
0.0931
0.2973
=pMm′
=pMm′
0.5001
0.1189
=pNt′
0.1588
0.1520
S11
p
(
N
y
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
N
t
′
=
p
M
m
+
O
x
′
;
p
N
t
′
=
p
N
t
+
L
i
′
10.5181
0.0024
11
1.6521
0.9478
0.0042
0.1006
0.3013
=pMm′
=pMm′
0.5203
0.1307
0.1590
=pNt′
0.1711
S12
p
(
Ny
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
N
t
′
=
p
M
m
+
O
x
′
;
p
N
t
′
=
p
N
t
+
M
m
′
9.2109
0.0045
11
1.5667
0.9461
0.0044
0.0927
0.2874
=pMm′
=pMm′
0.4388
0.1238
0.1494
0.1684
=pNt′
S13
p
(
N
y
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
N
t
′
=
p
M
m
+
O
x
′
;
p
N
t
′
<
p
N
t
+
P
j
′
0.0244
0.4473
10
1.7054
0.9780
0.0018
0.0928
0.2974
=pMm′
=pMm′
0.4972
0.1180
0.1750
=pNt+Ox′
=pNt+Ox′
S14
p
(
N
y
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0; pMm′=pMm+Pj′;pNt′<pNt+Pj′
0.0000
0.4528
9
1.7292
0.9755
0.0020
0.0929
0.3338
=pMm′
=pMm′
=pMm′
0.1171
0.1722
=pNt+Ox′
=pNt+Ox′
S15
p
(
Ny
)
(
C
,
i
,
B
,
a
i
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
P
j
′
;
p
N
t
′
=
p
N
t
+
P
j
′
7.6905
0.0097
8
1.5302
0.9224
0.0063
0.0981
0.3221
=pMm′
=pMm′
=pMm′
0.1244
=pNt′
=pNt′
=pNt′
S16
p
(
Ny
)
(
C
,
i
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
P
j
′
;
p
N
t
′
<
p
N
t
+
P
j
′
25.7016
0.0000
8
1.4843
0.0000
0.0151
0.0899
0.3239
=pMm′
=pMm′
=pMm′
0.0852
0.1631
=pNt+Ox′
=pNt+Ox′
S17
p
(
Ny
)
(
C
)
57.8675
0.0000
3
1.1455
0.0000
0.0813
0.2412
—
—
—
—
—
—
—
—
c
=
1.7180
The general model S1,p(Ny)(C,i,Bi,ai,M)all, specifies that: there is a constant unknown mortality of whitefly nymphs in all replicates (C); predation rates on whitefly nymphs varies with predator identity (i), blocks (Bi indicating that βi varies with predator identity) and whitefly nymph density (ai indicating that ai varies with predator identity); species-specific predation rates vary also between single- and multiple-predator treatment (M); 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→ pOx′=0); S2 → S3: no predation from P. phrygianus is further assumed (pOx′=0→pOx′=pLi′=0); S3 → S4: block effect constant between predator species (Bi → B); S4 → S5: effect of starting whitefly density on predation rates constant between predators (ai→a); S4 → S6: unknown mortality absent (C → ·); S4 → S7: per capita predation rates of M. melanotoma are identical in Mm and Mm + Nt treatments (pMm′=pMm+Nt′), and so on (S7 → S8 → S9; S8 → S10; S8 → S11; S8 → S12; S8 → S13 → S14 → S15; S14 → S16). Ki: number of model parameters (including the dispersion parameter of the negative binomial distribution, kNB, and overdispersion, ci). Model selection diagnostics are described in the method section. Parameter estimates are: unknown mortality (pCNy) and per capita predation rates of predators in single-(e.g., pMm′ for M. melanotoma) and multiple-predator treatments (e.g., pMm+Nt′ 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 countsEmerged, that is, the number of whitefly adults emerged in leaf samples. These observations were predicted as follows:
(7)Emerged(E)=sψ(B),
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 ψ(B) 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 countsNAd(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.Table 2
Best model structure and main parameter estimates relative to enemy biodiversity effects on whitefly nymph survival, adult emergence and survival (joint analysis).
Code
Model structure
Model selection diagnostics
Maximum likelihood estimates of model parameters
ΔQAICc
w
i
K
i
c
i
R
2
i
k
NB
Ψ
A
Ψ
B
p
O
x
′
p
O
x
+
N
t
′
p
O
x
+
L
i
′
p
O
x
+
M
m
′
p
L
i
′
p
L
i
+
N
t
′
p
L
i
+
O
x
′
p
L
i
+
M
m
′
J1
p
(
Ny
)(S13) Ψ(B) p(Ad)(C,i,Bi,ai,M)  all
64.6797
0.000
39
1.5351
0.1166
0.2077
0.4466
0.2272
0.2720
0.0000
0.1750
0.0000
0.5198
0.0000
0.0000
0.0000
J2
p
(
Ny
)(S13) Ψ(B) p(Ad)(i,Bi,ai,M)  all
61.5743
0.000
38
1.5255
0.1166
0.2077
0.4466
0.2272
0.2720
0.0000
0.1750
0.0000
0.5198
0.0000
0.0000
0.0000
J3
p
(
Ny
)(S13) Ψ(B) p(Ad)(i,M)  all
38.1612
0.000
30
1.4528
0.1166
0.2077
0.4473
0.2269
0.2725
0.0000
0.1751
0.0000
0.2745
0.0152
0.0000
0.0000
J4
p
(
Ny
)(S13) Ψ(B) p(Ad)(i,M);pNt′=0
27.2339
0.000
26
1.4190
0.1166
0.2077
0.4473
0.2269
0.2725
0.0000
0.0962
0.0000
0.2745
0.0815
0.0873
0.0000
J5
p
(
Ny
)(S13) Ψ(B) p(Ad)(i,M);pNt′=pMm′=0
16.8063
0.000
22
1.3867
0.1166
0.2077
0.4473
0.2269
0.2725
0.0000
0.0918
0.0000
0.2745
0.0815
0.0918
0.0000
J6
p
(
Ny
)(S13) Ψ(B) p(Ad)(i);pNt′=pMm′=0;pOx′=pOx+Pj′;pLi′=pLi+Pj′
5.6330
0.016
16
1.3440
0.0109
0.2326
0.4310
0.2234
0.0000
=pOx′
=pOx′
=pOx′
0.0498
=pLi′
=pLi′
=pLi′
J7
p
(
Ny
)(S13) Ψ(B) p(Ad)(i);pNt′=pMm′=0;pOx′=0;pLi′=pLi+Pj′
3.2648
0.051
15
1.3367
0.0109
0.2326
0.4310
0.2234
—
—
—
—
0.0498
=pLi′
=pLi′
=pLi′
J8
p
(
N
y
)
(S13)
Ψ
(
B
) p(Ad)(M);pNt'=pMm'=0;pOx′=pLi′
0.3907
0.216
15
1.3350
0.0886
0.2143
0.4411
0.2249
0.2642
—
—
—
=pOx′
—
—
—
J9
p
(
N
y
)
(S13)
Ψ
(
B
) p(Ad)(M);pNt'=pMm'=0;pOx'=pLi'=pLi+Ox'
0.0000
0.263
15
1.3351
0.1037
0.2108
0.4443
0.2266
0.2395
—
—
—
=pOx′
—
=pOx′
—
J10
p
(
N
y
)
(S13)
Ψ
(
B
) p(Ad)(M);pNt'=pMm'=0;  pOx'=pLi'=pOx+Li'=pLi+Ox'
1.0172
0.158
15
1.3362
0.0734
0.2179
0.4394
0.2258
0.1484
—
=pOx′
—
=pOx′
—
=pOx′
—
J11
p
(
N
y
)
(S13)
Ψ
(
B
) p(Ad)(k)
1.0771
0.153
14
1.3297
0.0000
0.2352
0.4275
0.2223
—
—
—
—
—
—
—
—
J12
p
(
Ny
)(S14) Ψ(B) p(Ad)(M);pNt′=pMm′=0;pOx′=pLi′=pLi+Ox′
1.2322
0.142
14
1.3455
0.1060
0.2102
0.4441
0.2264
0.2389
—
—
—
=pOx′
—
=pOx′
—
J13
p
(
Ny
)(S13) Ψ(B) p(Ad)(C)
3.4194
0.045
15
1.3370
0.0000
0.2352
0.4275
0.2223
—
—
—
—
—
—
—
—
J14
p
(
Ny
)(S13) Ψ(·) p(Ad)(M);pNt′=pMm′=0;pOx′=pLi′=pLi+Ox′
45.7700
0.000
14
1.3063
0.3911
0.3344
=ΨA
0.2036
—
—
—
=pOx′
—
=pOx′
—
c
=
1.3362
ML estimates of  parameters relative to predation on nymphs for Model J9
β
p
C
Ny
k
NB
a
N
t
a
M
m
p
M
m
′
p
M
m
+
N
t
′
p
M
m
+
O
x
′
p
M
m
+
L
i
′
p
N
t
′
p
N
t
+
O
x
′
p
N
t
+
L
i
′
p
N
t
+
M
m
′
−0.0892
0.0935
0.0019
0.0015
0.0076
0.2563
=pMm′
=pMm′
0.4286
0.1137
0.1686
=pNt+Ox′
=pNt+Ox′
The general model J1,p(Ny)(S13) Ψ(B) p(Ad)(C,i,Bi,ai,M)  all, specifies that: model structure for whitefly nymph survival is as in S13, p(Ny)(S13); emergence rates vary with blocks, Ψ(B); relatively to whitefly adult survival, p(Ad)(⋯), there is a constant unknown mortality in all replicates (C); predation varies with predator identity (i), blocks (Bi indicating that βi varies with predator identity) and whitefly density (ai indicating that ai varies with predator identity); species-specific predation rates vary also between single- and multiple-predator treatment (M); parameter structure refers to all four predators, (all). The variation of model structure following each simplification step is: J1 → J2: unknown mortality absent (C → ·); J2 → J3: block and starting whitefly density effects absent (Bi,ai → ·); J3 → J4: no predation from N. tenuis is assumed (all → pNt′=0); J4 → J5: no predation from M. melanotoma is further assumed (pNt′=0→pNt′=pMm′=0), and so on (J5 → J6; J5 → J7; J5 → J8; J5 → J9; J5 → J10; J9 → J12; J9 → J13). Ki: number of model parameters (including the dispersion parameter of the negative binomial distribution, kNB, and overdispersion, ci). Model selection diagnostics are described in the method section. Parameter estimates presented are: emergence rates for block A and B (ΨA and ΨB) and per capita predation rates of spiders in single- (e.g., pOx′ for O. lineatus) and multiple-predator treatments (e.g., pOx+Li′ for O. lineatuswhen 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 [31].The three series of counts (NNy(O) + Emerged + NAd(O); n=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.


#### 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 exponentsNi 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 Ni.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 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 Nt and Mm, 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 Nt and Mm, 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 Tables1 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 [51], was used to specify the error structure in the data [31, 33, 52]. With the negative binomial distribution, the variance = mean + kNB mean2. This distribution reduces to Poisson when kNB = 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 kNB 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 (c) of the corrected Akaike’s Information Criterion, QAICc, which is a modified version of the AIC used when the number of the estimable parameters in the model (Ki) is large relative to sample size (i.e., when n/K<40) and when data are overdispersed [30, 50];
(8)QAICc=-2log(ℒ(θ∣x,Mi))c+2Ki+2Ki(Ki+1)(n-Ki-1),
where n is the effective sample size, and log(ℒ(θ∣x,Mj)) is the log-likelihood evaluated at the maximum likelihood estimates, θ, given data x and model Mi. Overdispersion is estimated as follows:
(9)c^i=Deviancen-Ki,
where n is sample size and Ki is the number of estimated parameters of model Mi [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≥c≤4 [30]. Overdispersion (ci), as well as the negative binomial kNB, is estimated as a parameter and is therefore always considered in the Ki, the number of estimated parameters of model Mi. Moreover, the QAICc 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], c=∑wic^i.Finally, support for each modelMi was assessed by using ΔQAICc=QAICc(i)-LowestQAICc and the Akaike weight, wi, and a measure of the explanatory power (R2) of the best models selected was estimated using the dispersion parameter kNB of the following negative binomial [53]:
(10)R2=1-kbestkmax,
where kmax is estimated in the negative binomial model with only a constant mortality term, p(C).


## 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 whiteflyTrialeurodes 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–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].


## 2.2. Experimental Design

The effect of single-predator (N. tenuis, M. melanotoma, O. lineatus,and P. phrygianus;thereafter, Nt, Mm, Ox, andLi, resp.) and two-predator (Nt + Mm, Nt + Ox, Nt + Li, Mm + Ox, Mm + Li,andOx + 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 [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.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 inspection. Insects were stored in alcohol and subsequently counted by means of a stereomicroscope. The number of whitefly adults (NAd(O)) and predators (NNt, NMm, NLi, and NOx) 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 (NNy-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 (NNy(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.


## 2.3. Modeling Enemy Biodiversity Effects on Prey

The starting point of model selection approaches is to develop a set ofa 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 (NNy) and adults (NAd) present in each experimental unit could be approximated by the following:
(1)NNy≈NNy-EU,NAd≈0.During the experiment, the relative numbers of the two whitefly subpopulations change as a consequence of emergence rates, ψ(B), and predation. Predation results in a reduction of the two subpopulations proportional to predator density, Ni, and predator-specific per capita predation rates on whitefly nymphs, pi(Ny)′, and adults, pi(Ad)′, where i is the predator species.If the number of whiteflies surviving predation from a given predator speciesi is approximated by the following general multiplicative risk model [25, 26, 45–47]:
(2)NWhitefly(t+1)=NWhitefly(t)(1-pi′)Ni,
the number of whitefly adults expected at the end of the experiment can be written as follows:
(3)NAd(E)=NNy-EU(1-pNt(Ny)′)NNt(1-pMm(Ny)′)NMm×(1-pOx(Ny)′)NOx(1-pLi(Ny)′)NLiψ(B)×(1-pNt(Ad)′)NNt(1-pMm(Ad)′)Nm×(1-pOx(Ad)′)NOx(1-pLi(Ad)′)NLi.
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 (NNy-EU) that survive predation on the nymph stage due to different predators: NNy-EU∏in‍(1-pi(Ny)′)Ni, where i is predator species Nt, Mm, Ox, or Li. Surviving nymphs become adults with rate ψ(B), and to be counted at the end of the experiment, they have to survive predation on the adult stage (∏in‍(1-pi(Ad)′)Ni). 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 Table5) 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 (NNy-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; n=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:
(4)pi(block,unit)′=pi′+βie(NNy-EU*ai),
where βi is an additive factor that can vary with predator identity and applies only to block A, whereas ai is a scaling factor that can vary with predator identity and relates per capita predation rate with whitefly nymph density NNy-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:(5)pi(Ny)′≠pi+Pj(Ny)′,pi(Ad)′≠pi+Pj(Ad)′,
where Pj refers to the presence of another predator j. 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 (pCAd) and nymph subpopulations (pCNy) 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) [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 (NNy(O)), whitefly adult emergence (Emerged), and whitefly adult survival (NAd(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.


### 2.4.1. Single Series Analysis

The first part of the process—whitefly nymph survival—was studied by analysing the series of countsNNy(O), that is, the number of whitefly nymphs survived in leaf samples. These observations were predicted as follows:
(6)NNy(E)=s(1-pNt(Ny)′)NNt(1-pMm(Ny)′)NMm×(1-pOx(Ny)′)NOx(1-pLi(Ny)′)NLi(1-pCNy),
where NNy(E) is the expected number of whitefly nymphs present in leaf samples that survived at the end of the experiment and s 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 (+kNB and c = 27, see below), was simplified according to the details presented in Table 1.Table 1
Best model structure and main parameter estimates relative to enemy biodiversity effects on whitefly nymph survival.
Code
Model structure
Model selection diagnostics
Maximum likelihood estimates of model parameters
ΔQAICc
w
i
K
i
c
i
R
2
i
k
NB
p
C
Ny
p
M
m
′
p
M
m
+
N
t
′
p
M
m
+
O
x
′
p
M
m
+
L
i
′
p
N
t
′
p
N
t
+
O
x
′
p
N
t
+
L
i
′
p
N
t
+
M
m
′
S1
p
(
Ny
)
(
C
,
i
,
B
i
,
a
i
,
M
);  all
57.6253
0.0000
27
2.3675
1.0000
0.0000
0.0729
0.3004
0.0000
0.0000
0.4914
0.1611
0.0256
0.1192
0.2908
S2
p
(
Ny
)
(
C
,
i
,
B
i
,
a
i
,
M
)
;
p
O
x
′
=
0
34.4510
0.0000
21
2.2286
0.9944
0.0005
0.0869
0.2860
0.4184
0.2474
0.4899
0.1097
0.1541
0.1787
0.1552
S3
p
(
Ny
)
(
C
,
i
,
B
i
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
14.2799
0.0004
15
1.8687
0.9810
0.0015
0.0894
0.3128
0.3905
0.2870
0.5029
0.1163
0.1633
0.1802
0.1700
S4
p
(
Ny
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
11.5475
0.0014
14
1.8225
0.9787
0.0017
0.0914
0.3001
0.3798
0.2615
0.4784
0.1148
0.1614
0.1785
0.1678
S5
p
(
Ny
)
(
C
,
i
,
B
,
a
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
26.3280
0.0000
13
1.6959
0.9072
0.0075
0.0906
0.1021
0.1645
0.0850
0.1304
0.1336
0.2037
0.2306
0.1646
S6
p
(
Ny
)
(
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
24.1824
0.0000
13
1.4922
0.8872
0.0092
—
0.2716
0.2591
0.2200
0.3759
0.1208
0.1598
0.1878
0.1662
S7
p
(
Ny
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
N
t
′
8.3366
0.0070
13
1.7840
0.9783
0.0018
0.0914
0.3019
=pMm′
0.2616
0.4788
0.1159
0.1633
0.1804
0.1739
S8
p
(
Ny
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
N
t
′
=
p
M
m
+
O
x
′
5.3157
0.0317
12
1.7565
0.9785
0.0017
0.0917
0.2976
=pMm′
=pMm′
0.5006
0.1162
0.1638
0.1809
0.1757
S9
p
(
Ny
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
P
j
′
5.1225
0.0350
11
1.7795
0.9759
0.0020
0.0917
0.3323
=pMm′
=pMm′
=pMm′
0.1152
0.1619
0.1792
0.1704
S10
p
(
Ny
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
N
t
′
=
p
M
m
+
O
x
′
;
p
N
t
′
=
p
N
t
+
O
x
′
8.1180
0.0078
11
1.6109
0.9515
0.0039
0.0931
0.2973
=pMm′
=pMm′
0.5001
0.1189
=pNt′
0.1588
0.1520
S11
p
(
N
y
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
N
t
′
=
p
M
m
+
O
x
′
;
p
N
t
′
=
p
N
t
+
L
i
′
10.5181
0.0024
11
1.6521
0.9478
0.0042
0.1006
0.3013
=pMm′
=pMm′
0.5203
0.1307
0.1590
=pNt′
0.1711
S12
p
(
Ny
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
N
t
′
=
p
M
m
+
O
x
′
;
p
N
t
′
=
p
N
t
+
M
m
′
9.2109
0.0045
11
1.5667
0.9461
0.0044
0.0927
0.2874
=pMm′
=pMm′
0.4388
0.1238
0.1494
0.1684
=pNt′
S13
p
(
N
y
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
N
t
′
=
p
M
m
+
O
x
′
;
p
N
t
′
<
p
N
t
+
P
j
′
0.0244
0.4473
10
1.7054
0.9780
0.0018
0.0928
0.2974
=pMm′
=pMm′
0.4972
0.1180
0.1750
=pNt+Ox′
=pNt+Ox′
S14
p
(
N
y
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0; pMm′=pMm+Pj′;pNt′<pNt+Pj′
0.0000
0.4528
9
1.7292
0.9755
0.0020
0.0929
0.3338
=pMm′
=pMm′
=pMm′
0.1171
0.1722
=pNt+Ox′
=pNt+Ox′
S15
p
(
Ny
)
(
C
,
i
,
B
,
a
i
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
P
j
′
;
p
N
t
′
=
p
N
t
+
P
j
′
7.6905
0.0097
8
1.5302
0.9224
0.0063
0.0981
0.3221
=pMm′
=pMm′
=pMm′
0.1244
=pNt′
=pNt′
=pNt′
S16
p
(
Ny
)
(
C
,
i
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
P
j
′
;
p
N
t
′
<
p
N
t
+
P
j
′
25.7016
0.0000
8
1.4843
0.0000
0.0151
0.0899
0.3239
=pMm′
=pMm′
=pMm′
0.0852
0.1631
=pNt+Ox′
=pNt+Ox′
S17
p
(
Ny
)
(
C
)
57.8675
0.0000
3
1.1455
0.0000
0.0813
0.2412
—
—
—
—
—
—
—
—
c
=
1.7180
The general model S1,p(Ny)(C,i,Bi,ai,M)all, specifies that: there is a constant unknown mortality of whitefly nymphs in all replicates (C); predation rates on whitefly nymphs varies with predator identity (i), blocks (Bi indicating that βi varies with predator identity) and whitefly nymph density (ai indicating that ai varies with predator identity); species-specific predation rates vary also between single- and multiple-predator treatment (M); 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→ pOx′=0); S2 → S3: no predation from P. phrygianus is further assumed (pOx′=0→pOx′=pLi′=0); S3 → S4: block effect constant between predator species (Bi → B); S4 → S5: effect of starting whitefly density on predation rates constant between predators (ai→a); S4 → S6: unknown mortality absent (C → ·); S4 → S7: per capita predation rates of M. melanotoma are identical in Mm and Mm + Nt treatments (pMm′=pMm+Nt′), and so on (S7 → S8 → S9; S8 → S10; S8 → S11; S8 → S12; S8 → S13 → S14 → S15; S14 → S16). Ki: number of model parameters (including the dispersion parameter of the negative binomial distribution, kNB, and overdispersion, ci). Model selection diagnostics are described in the method section. Parameter estimates are: unknown mortality (pCNy) and per capita predation rates of predators in single-(e.g., pMm′ for M. melanotoma) and multiple-predator treatments (e.g., pMm+Nt′ 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 countsEmerged, that is, the number of whitefly adults emerged in leaf samples. These observations were predicted as follows:
(7)Emerged(E)=sψ(B),
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 ψ(B) 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 countsNAd(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.Table 2
Best model structure and main parameter estimates relative to enemy biodiversity effects on whitefly nymph survival, adult emergence and survival (joint analysis).
Code
Model structure
Model selection diagnostics
Maximum likelihood estimates of model parameters
ΔQAICc
w
i
K
i
c
i
R
2
i
k
NB
Ψ
A
Ψ
B
p
O
x
′
p
O
x
+
N
t
′
p
O
x
+
L
i
′
p
O
x
+
M
m
′
p
L
i
′
p
L
i
+
N
t
′
p
L
i
+
O
x
′
p
L
i
+
M
m
′
J1
p
(
Ny
)(S13) Ψ(B) p(Ad)(C,i,Bi,ai,M)  all
64.6797
0.000
39
1.5351
0.1166
0.2077
0.4466
0.2272
0.2720
0.0000
0.1750
0.0000
0.5198
0.0000
0.0000
0.0000
J2
p
(
Ny
)(S13) Ψ(B) p(Ad)(i,Bi,ai,M)  all
61.5743
0.000
38
1.5255
0.1166
0.2077
0.4466
0.2272
0.2720
0.0000
0.1750
0.0000
0.5198
0.0000
0.0000
0.0000
J3
p
(
Ny
)(S13) Ψ(B) p(Ad)(i,M)  all
38.1612
0.000
30
1.4528
0.1166
0.2077
0.4473
0.2269
0.2725
0.0000
0.1751
0.0000
0.2745
0.0152
0.0000
0.0000
J4
p
(
Ny
)(S13) Ψ(B) p(Ad)(i,M);pNt′=0
27.2339
0.000
26
1.4190
0.1166
0.2077
0.4473
0.2269
0.2725
0.0000
0.0962
0.0000
0.2745
0.0815
0.0873
0.0000
J5
p
(
Ny
)(S13) Ψ(B) p(Ad)(i,M);pNt′=pMm′=0
16.8063
0.000
22
1.3867
0.1166
0.2077
0.4473
0.2269
0.2725
0.0000
0.0918
0.0000
0.2745
0.0815
0.0918
0.0000
J6
p
(
Ny
)(S13) Ψ(B) p(Ad)(i);pNt′=pMm′=0;pOx′=pOx+Pj′;pLi′=pLi+Pj′
5.6330
0.016
16
1.3440
0.0109
0.2326
0.4310
0.2234
0.0000
=pOx′
=pOx′
=pOx′
0.0498
=pLi′
=pLi′
=pLi′
J7
p
(
Ny
)(S13) Ψ(B) p(Ad)(i);pNt′=pMm′=0;pOx′=0;pLi′=pLi+Pj′
3.2648
0.051
15
1.3367
0.0109
0.2326
0.4310
0.2234
—
—
—
—
0.0498
=pLi′
=pLi′
=pLi′
J8
p
(
N
y
)
(S13)
Ψ
(
B
) p(Ad)(M);pNt'=pMm'=0;pOx′=pLi′
0.3907
0.216
15
1.3350
0.0886
0.2143
0.4411
0.2249
0.2642
—
—
—
=pOx′
—
—
—
J9
p
(
N
y
)
(S13)
Ψ
(
B
) p(Ad)(M);pNt'=pMm'=0;pOx'=pLi'=pLi+Ox'
0.0000
0.263
15
1.3351
0.1037
0.2108
0.4443
0.2266
0.2395
—
—
—
=pOx′
—
=pOx′
—
J10
p
(
N
y
)
(S13)
Ψ
(
B
) p(Ad)(M);pNt'=pMm'=0;  pOx'=pLi'=pOx+Li'=pLi+Ox'
1.0172
0.158
15
1.3362
0.0734
0.2179
0.4394
0.2258
0.1484
—
=pOx′
—
=pOx′
—
=pOx′
—
J11
p
(
N
y
)
(S13)
Ψ
(
B
) p(Ad)(k)
1.0771
0.153
14
1.3297
0.0000
0.2352
0.4275
0.2223
—
—
—
—
—
—
—
—
J12
p
(
Ny
)(S14) Ψ(B) p(Ad)(M);pNt′=pMm′=0;pOx′=pLi′=pLi+Ox′
1.2322
0.142
14
1.3455
0.1060
0.2102
0.4441
0.2264
0.2389
—
—
—
=pOx′
—
=pOx′
—
J13
p
(
Ny
)(S13) Ψ(B) p(Ad)(C)
3.4194
0.045
15
1.3370
0.0000
0.2352
0.4275
0.2223
—
—
—
—
—
—
—
—
J14
p
(
Ny
)(S13) Ψ(·) p(Ad)(M);pNt′=pMm′=0;pOx′=pLi′=pLi+Ox′
45.7700
0.000
14
1.3063
0.3911
0.3344
=ΨA
0.2036
—
—
—
=pOx′
—
=pOx′
—
c
=
1.3362
ML estimates of  parameters relative to predation on nymphs for Model J9
β
p
C
Ny
k
NB
a
N
t
a
M
m
p
M
m
′
p
M
m
+
N
t
′
p
M
m
+
O
x
′
p
M
m
+
L
i
′
p
N
t
′
p
N
t
+
O
x
′
p
N
t
+
L
i
′
p
N
t
+
M
m
′
−0.0892
0.0935
0.0019
0.0015
0.0076
0.2563
=pMm′
=pMm′
0.4286
0.1137
0.1686
=pNt+Ox′
=pNt+Ox′
The general model J1,p(Ny)(S13) Ψ(B) p(Ad)(C,i,Bi,ai,M)  all, specifies that: model structure for whitefly nymph survival is as in S13, p(Ny)(S13); emergence rates vary with blocks, Ψ(B); relatively to whitefly adult survival, p(Ad)(⋯), there is a constant unknown mortality in all replicates (C); predation varies with predator identity (i), blocks (Bi indicating that βi varies with predator identity) and whitefly density (ai indicating that ai varies with predator identity); species-specific predation rates vary also between single- and multiple-predator treatment (M); parameter structure refers to all four predators, (all). The variation of model structure following each simplification step is: J1 → J2: unknown mortality absent (C → ·); J2 → J3: block and starting whitefly density effects absent (Bi,ai → ·); J3 → J4: no predation from N. tenuis is assumed (all → pNt′=0); J4 → J5: no predation from M. melanotoma is further assumed (pNt′=0→pNt′=pMm′=0), and so on (J5 → J6; J5 → J7; J5 → J8; J5 → J9; J5 → J10; J9 → J12; J9 → J13). Ki: number of model parameters (including the dispersion parameter of the negative binomial distribution, kNB, and overdispersion, ci). Model selection diagnostics are described in the method section. Parameter estimates presented are: emergence rates for block A and B (ΨA and ΨB) and per capita predation rates of spiders in single- (e.g., pOx′ for O. lineatus) and multiple-predator treatments (e.g., pOx+Li′ for O. lineatuswhen 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 [31].The three series of counts (NNy(O) + Emerged + NAd(O); n=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.


### 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 exponentsNi 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 Ni.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 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 Nt and Mm, 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 Nt and Mm, 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 Tables1 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 [51], was used to specify the error structure in the data [31, 33, 52]. With the negative binomial distribution, the variance = mean + kNB mean2. This distribution reduces to Poisson when kNB = 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 kNB 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 (c) of the corrected Akaike’s Information Criterion, QAICc, which is a modified version of the AIC used when the number of the estimable parameters in the model (Ki) is large relative to sample size (i.e., when n/K<40) and when data are overdispersed [30, 50];
(8)QAICc=-2log(ℒ(θ∣x,Mi))c+2Ki+2Ki(Ki+1)(n-Ki-1),
where n is the effective sample size, and log(ℒ(θ∣x,Mj)) is the log-likelihood evaluated at the maximum likelihood estimates, θ, given data x and model Mi. Overdispersion is estimated as follows:
(9)c^i=Deviancen-Ki,
where n is sample size and Ki is the number of estimated parameters of model Mi [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≥c≤4 [30]. Overdispersion (ci), as well as the negative binomial kNB, is estimated as a parameter and is therefore always considered in the Ki, the number of estimated parameters of model Mi. Moreover, the QAICc 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], c=∑wic^i.Finally, support for each modelMi was assessed by using ΔQAICc=QAICc(i)-LowestQAICc and the Akaike weight, wi, and a measure of the explanatory power (R2) of the best models selected was estimated using the dispersion parameter kNB of the following negative binomial [53]:
(10)R2=1-kbestkmax,
where kmax is estimated in the negative binomial model with only a constant mortality term, p(C).


## 2.4.1. Single Series Analysis

The first part of the process—whitefly nymph survival—was studied by analysing the series of countsNNy(O), that is, the number of whitefly nymphs survived in leaf samples. These observations were predicted as follows:
(6)NNy(E)=s(1-pNt(Ny)′)NNt(1-pMm(Ny)′)NMm×(1-pOx(Ny)′)NOx(1-pLi(Ny)′)NLi(1-pCNy),
where NNy(E) is the expected number of whitefly nymphs present in leaf samples that survived at the end of the experiment and s 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 (+kNB and c = 27, see below), was simplified according to the details presented in Table 1.Table 1
Best model structure and main parameter estimates relative to enemy biodiversity effects on whitefly nymph survival.
Code
Model structure
Model selection diagnostics
Maximum likelihood estimates of model parameters
ΔQAICc
w
i
K
i
c
i
R
2
i
k
NB
p
C
Ny
p
M
m
′
p
M
m
+
N
t
′
p
M
m
+
O
x
′
p
M
m
+
L
i
′
p
N
t
′
p
N
t
+
O
x
′
p
N
t
+
L
i
′
p
N
t
+
M
m
′
S1
p
(
Ny
)
(
C
,
i
,
B
i
,
a
i
,
M
);  all
57.6253
0.0000
27
2.3675
1.0000
0.0000
0.0729
0.3004
0.0000
0.0000
0.4914
0.1611
0.0256
0.1192
0.2908
S2
p
(
Ny
)
(
C
,
i
,
B
i
,
a
i
,
M
)
;
p
O
x
′
=
0
34.4510
0.0000
21
2.2286
0.9944
0.0005
0.0869
0.2860
0.4184
0.2474
0.4899
0.1097
0.1541
0.1787
0.1552
S3
p
(
Ny
)
(
C
,
i
,
B
i
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
14.2799
0.0004
15
1.8687
0.9810
0.0015
0.0894
0.3128
0.3905
0.2870
0.5029
0.1163
0.1633
0.1802
0.1700
S4
p
(
Ny
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
11.5475
0.0014
14
1.8225
0.9787
0.0017
0.0914
0.3001
0.3798
0.2615
0.4784
0.1148
0.1614
0.1785
0.1678
S5
p
(
Ny
)
(
C
,
i
,
B
,
a
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
26.3280
0.0000
13
1.6959
0.9072
0.0075
0.0906
0.1021
0.1645
0.0850
0.1304
0.1336
0.2037
0.2306
0.1646
S6
p
(
Ny
)
(
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
24.1824
0.0000
13
1.4922
0.8872
0.0092
—
0.2716
0.2591
0.2200
0.3759
0.1208
0.1598
0.1878
0.1662
S7
p
(
Ny
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
N
t
′
8.3366
0.0070
13
1.7840
0.9783
0.0018
0.0914
0.3019
=pMm′
0.2616
0.4788
0.1159
0.1633
0.1804
0.1739
S8
p
(
Ny
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
N
t
′
=
p
M
m
+
O
x
′
5.3157
0.0317
12
1.7565
0.9785
0.0017
0.0917
0.2976
=pMm′
=pMm′
0.5006
0.1162
0.1638
0.1809
0.1757
S9
p
(
Ny
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
P
j
′
5.1225
0.0350
11
1.7795
0.9759
0.0020
0.0917
0.3323
=pMm′
=pMm′
=pMm′
0.1152
0.1619
0.1792
0.1704
S10
p
(
Ny
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
N
t
′
=
p
M
m
+
O
x
′
;
p
N
t
′
=
p
N
t
+
O
x
′
8.1180
0.0078
11
1.6109
0.9515
0.0039
0.0931
0.2973
=pMm′
=pMm′
0.5001
0.1189
=pNt′
0.1588
0.1520
S11
p
(
N
y
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
N
t
′
=
p
M
m
+
O
x
′
;
p
N
t
′
=
p
N
t
+
L
i
′
10.5181
0.0024
11
1.6521
0.9478
0.0042
0.1006
0.3013
=pMm′
=pMm′
0.5203
0.1307
0.1590
=pNt′
0.1711
S12
p
(
Ny
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
N
t
′
=
p
M
m
+
O
x
′
;
p
N
t
′
=
p
N
t
+
M
m
′
9.2109
0.0045
11
1.5667
0.9461
0.0044
0.0927
0.2874
=pMm′
=pMm′
0.4388
0.1238
0.1494
0.1684
=pNt′
S13
p
(
N
y
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
N
t
′
=
p
M
m
+
O
x
′
;
p
N
t
′
<
p
N
t
+
P
j
′
0.0244
0.4473
10
1.7054
0.9780
0.0018
0.0928
0.2974
=pMm′
=pMm′
0.4972
0.1180
0.1750
=pNt+Ox′
=pNt+Ox′
S14
p
(
N
y
)
(
C
,
i
,
B
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0; pMm′=pMm+Pj′;pNt′<pNt+Pj′
0.0000
0.4528
9
1.7292
0.9755
0.0020
0.0929
0.3338
=pMm′
=pMm′
=pMm′
0.1171
0.1722
=pNt+Ox′
=pNt+Ox′
S15
p
(
Ny
)
(
C
,
i
,
B
,
a
i
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
P
j
′
;
p
N
t
′
=
p
N
t
+
P
j
′
7.6905
0.0097
8
1.5302
0.9224
0.0063
0.0981
0.3221
=pMm′
=pMm′
=pMm′
0.1244
=pNt′
=pNt′
=pNt′
S16
p
(
Ny
)
(
C
,
i
,
a
i
,
M
)
;
p
O
x
′
=
p
L
i
′
=
0
;
p
M
m
′
=
p
M
m
+
P
j
′
;
p
N
t
′
<
p
N
t
+
P
j
′
25.7016
0.0000
8
1.4843
0.0000
0.0151
0.0899
0.3239
=pMm′
=pMm′
=pMm′
0.0852
0.1631
=pNt+Ox′
=pNt+Ox′
S17
p
(
Ny
)
(
C
)
57.8675
0.0000
3
1.1455
0.0000
0.0813
0.2412
—
—
—
—
—
—
—
—
c
=
1.7180
The general model S1,p(Ny)(C,i,Bi,ai,M)all, specifies that: there is a constant unknown mortality of whitefly nymphs in all replicates (C); predation rates on whitefly nymphs varies with predator identity (i), blocks (Bi indicating that βi varies with predator identity) and whitefly nymph density (ai indicating that ai varies with predator identity); species-specific predation rates vary also between single- and multiple-predator treatment (M); 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→ pOx′=0); S2 → S3: no predation from P. phrygianus is further assumed (pOx′=0→pOx′=pLi′=0); S3 → S4: block effect constant between predator species (Bi → B); S4 → S5: effect of starting whitefly density on predation rates constant between predators (ai→a); S4 → S6: unknown mortality absent (C → ·); S4 → S7: per capita predation rates of M. melanotoma are identical in Mm and Mm + Nt treatments (pMm′=pMm+Nt′), and so on (S7 → S8 → S9; S8 → S10; S8 → S11; S8 → S12; S8 → S13 → S14 → S15; S14 → S16). Ki: number of model parameters (including the dispersion parameter of the negative binomial distribution, kNB, and overdispersion, ci). Model selection diagnostics are described in the method section. Parameter estimates are: unknown mortality (pCNy) and per capita predation rates of predators in single-(e.g., pMm′ for M. melanotoma) and multiple-predator treatments (e.g., pMm+Nt′ 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 countsEmerged, that is, the number of whitefly adults emerged in leaf samples. These observations were predicted as follows:
(7)Emerged(E)=sψ(B),
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 ψ(B) 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 countsNAd(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.Table 2
Best model structure and main parameter estimates relative to enemy biodiversity effects on whitefly nymph survival, adult emergence and survival (joint analysis).
Code
Model structure
Model selection diagnostics
Maximum likelihood estimates of model parameters
ΔQAICc
w
i
K
i
c
i
R
2
i
k
NB
Ψ
A
Ψ
B
p
O
x
′
p
O
x
+
N
t
′
p
O
x
+
L
i
′
p
O
x
+
M
m
′
p
L
i
′
p
L
i
+
N
t
′
p
L
i
+
O
x
′
p
L
i
+
M
m
′
J1
p
(
Ny
)(S13) Ψ(B) p(Ad)(C,i,Bi,ai,M)  all
64.6797
0.000
39
1.5351
0.1166
0.2077
0.4466
0.2272
0.2720
0.0000
0.1750
0.0000
0.5198
0.0000
0.0000
0.0000
J2
p
(
Ny
)(S13) Ψ(B) p(Ad)(i,Bi,ai,M)  all
61.5743
0.000
38
1.5255
0.1166
0.2077
0.4466
0.2272
0.2720
0.0000
0.1750
0.0000
0.5198
0.0000
0.0000
0.0000
J3
p
(
Ny
)(S13) Ψ(B) p(Ad)(i,M)  all
38.1612
0.000
30
1.4528
0.1166
0.2077
0.4473
0.2269
0.2725
0.0000
0.1751
0.0000
0.2745
0.0152
0.0000
0.0000
J4
p
(
Ny
)(S13) Ψ(B) p(Ad)(i,M);pNt′=0
27.2339
0.000
26
1.4190
0.1166
0.2077
0.4473
0.2269
0.2725
0.0000
0.0962
0.0000
0.2745
0.0815
0.0873
0.0000
J5
p
(
Ny
)(S13) Ψ(B) p(Ad)(i,M);pNt′=pMm′=0
16.8063
0.000
22
1.3867
0.1166
0.2077
0.4473
0.2269
0.2725
0.0000
0.0918
0.0000
0.2745
0.0815
0.0918
0.0000
J6
p
(
Ny
)(S13) Ψ(B) p(Ad)(i);pNt′=pMm′=0;pOx′=pOx+Pj′;pLi′=pLi+Pj′
5.6330
0.016
16
1.3440
0.0109
0.2326
0.4310
0.2234
0.0000
=pOx′
=pOx′
=pOx′
0.0498
=pLi′
=pLi′
=pLi′
J7
p
(
Ny
)(S13) Ψ(B) p(Ad)(i);pNt′=pMm′=0;pOx′=0;pLi′=pLi+Pj′
3.2648
0.051
15
1.3367
0.0109
0.2326
0.4310
0.2234
—
—
—
—
0.0498
=pLi′
=pLi′
=pLi′
J8
p
(
N
y
)
(S13)
Ψ
(
B
) p(Ad)(M);pNt'=pMm'=0;pOx′=pLi′
0.3907
0.216
15
1.3350
0.0886
0.2143
0.4411
0.2249
0.2642
—
—
—
=pOx′
—
—
—
J9
p
(
N
y
)
(S13)
Ψ
(
B
) p(Ad)(M);pNt'=pMm'=0;pOx'=pLi'=pLi+Ox'
0.0000
0.263
15
1.3351
0.1037
0.2108
0.4443
0.2266
0.2395
—
—
—
=pOx′
—
=pOx′
—
J10
p
(
N
y
)
(S13)
Ψ
(
B
) p(Ad)(M);pNt'=pMm'=0;  pOx'=pLi'=pOx+Li'=pLi+Ox'
1.0172
0.158
15
1.3362
0.0734
0.2179
0.4394
0.2258
0.1484
—
=pOx′
—
=pOx′
—
=pOx′
—
J11
p
(
N
y
)
(S13)
Ψ
(
B
) p(Ad)(k)
1.0771
0.153
14
1.3297
0.0000
0.2352
0.4275
0.2223
—
—
—
—
—
—
—
—
J12
p
(
Ny
)(S14) Ψ(B) p(Ad)(M);pNt′=pMm′=0;pOx′=pLi′=pLi+Ox′
1.2322
0.142
14
1.3455
0.1060
0.2102
0.4441
0.2264
0.2389
—
—
—
=pOx′
—
=pOx′
—
J13
p
(
Ny
)(S13) Ψ(B) p(Ad)(C)
3.4194
0.045
15
1.3370
0.0000
0.2352
0.4275
0.2223
—
—
—
—
—
—
—
—
J14
p
(
Ny
)(S13) Ψ(·) p(Ad)(M);pNt′=pMm′=0;pOx′=pLi′=pLi+Ox′
45.7700
0.000
14
1.3063
0.3911
0.3344
=ΨA
0.2036
—
—
—
=pOx′
—
=pOx′
—
c
=
1.3362
ML estimates of  parameters relative to predation on nymphs for Model J9
β
p
C
Ny
k
NB
a
N
t
a
M
m
p
M
m
′
p
M
m
+
N
t
′
p
M
m
+
O
x
′
p
M
m
+
L
i
′
p
N
t
′
p
N
t
+
O
x
′
p
N
t
+
L
i
′
p
N
t
+
M
m
′
−0.0892
0.0935
0.0019
0.0015
0.0076
0.2563
=pMm′
=pMm′
0.4286
0.1137
0.1686
=pNt+Ox′
=pNt+Ox′
The general model J1,p(Ny)(S13) Ψ(B) p(Ad)(C,i,Bi,ai,M)  all, specifies that: model structure for whitefly nymph survival is as in S13, p(Ny)(S13); emergence rates vary with blocks, Ψ(B); relatively to whitefly adult survival, p(Ad)(⋯), there is a constant unknown mortality in all replicates (C); predation varies with predator identity (i), blocks (Bi indicating that βi varies with predator identity) and whitefly density (ai indicating that ai varies with predator identity); species-specific predation rates vary also between single- and multiple-predator treatment (M); parameter structure refers to all four predators, (all). The variation of model structure following each simplification step is: J1 → J2: unknown mortality absent (C → ·); J2 → J3: block and starting whitefly density effects absent (Bi,ai → ·); J3 → J4: no predation from N. tenuis is assumed (all → pNt′=0); J4 → J5: no predation from M. melanotoma is further assumed (pNt′=0→pNt′=pMm′=0), and so on (J5 → J6; J5 → J7; J5 → J8; J5 → J9; J5 → J10; J9 → J12; J9 → J13). Ki: number of model parameters (including the dispersion parameter of the negative binomial distribution, kNB, and overdispersion, ci). Model selection diagnostics are described in the method section. Parameter estimates presented are: emergence rates for block A and B (ΨA and ΨB) and per capita predation rates of spiders in single- (e.g., pOx′ for O. lineatus) and multiple-predator treatments (e.g., pOx+Li′ for O. lineatuswhen 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 [31].The three series of counts (NNy(O) + Emerged + NAd(O); n=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.


## 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 exponentsNi 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 Ni.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 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 Nt and Mm, 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 Nt and Mm, 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 Tables1 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 [51], was used to specify the error structure in the data [31, 33, 52]. With the negative binomial distribution, the variance = mean + kNB mean2. This distribution reduces to Poisson when kNB = 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 kNB 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 (c) of the corrected Akaike’s Information Criterion, QAICc, which is a modified version of the AIC used when the number of the estimable parameters in the model (Ki) is large relative to sample size (i.e., when n/K<40) and when data are overdispersed [30, 50];
(8)QAICc=-2log(ℒ(θ∣x,Mi))c+2Ki+2Ki(Ki+1)(n-Ki-1),
where n is the effective sample size, and log(ℒ(θ∣x,Mj)) is the log-likelihood evaluated at the maximum likelihood estimates, θ, given data x and model Mi. Overdispersion is estimated as follows:
(9)c^i=Deviancen-Ki,
where n is sample size and Ki is the number of estimated parameters of model Mi [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≥c≤4 [30]. Overdispersion (ci), as well as the negative binomial kNB, is estimated as a parameter and is therefore always considered in the Ki, the number of estimated parameters of model Mi. Moreover, the QAICc 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], c=∑wic^i.Finally, support for each modelMi was assessed by using ΔQAICc=QAICc(i)-LowestQAICc and the Akaike weight, wi, and a measure of the explanatory power (R2) of the best models selected was estimated using the dispersion parameter kNB of the following negative binomial [53]:
(10)R2=1-kbestkmax,
where kmax is estimated in the negative binomial model with only a constant mortality term, p(C).


## 3. Results

### 3.1. What Do Raw Data Suggest?

Table5 shows how counts of the 3 data series analysed varied among experimental units and between blocks: NNy(O) = 238±81 (SD) and 121±69 for blocks A and B, respectively; Emerged = 110±35 and 39±17 for A and B; NAd(O) = 263±142 and 60±42 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 (NNy-EU = 655±323 and 427±262; 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 NAd(O)/NNy-EU = 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 Table3 (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.Table 3
Mean mortality of whitefly nymphs (Dead/s) for each treatment and control (n=6).
Nt
Mm
Ox
Li
Nt
0.4187
Mm
0.4016
0.2562
Ox
0.4009
0.2009
0.0673
Li
0.4159
0.2438
0.1042
0.0792
C
0.0978
Dead: number of dead whitefly nymphs counted in leaf samples; s: sample size (Dead + Alive + Emerged); Nt, Mm, Ox, Li: predator treatments (N. tenuis, M. melanotoma, O. lineatus, and P. phrygianus resp.), C: control.Table 4
Mean number of predators recovered in each treatment (±SD;n=6).
Nt
Mm
Ox
Li
Nt
8.83
±
5.95
3.50
±
2.07
1.33
±
0.82
1.33
±
1.51
Mm
5.50
±
4.04
6.17
±
1.33
1.33
±
0.82
1.50
±
1.05
Ox
8.00
±
4.34
7.17
±
2.79
0.67
±
0.52
1.17
±
1.17
Li
5.67
±
3.83
5.67
±
2.50
0.50
±
0.55
1.67
±
0.82
Mean
7.00
5.63
0.96
1.42
Number of predators recovered for each experimental unit are given in Table5 (NNt, NMm, NOx, NLi), 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.Table 5
Number of whitefly nymphs, adults and predators counted at the end of the experiment and model expectations.
Treatment
Block
Counts in Experimental Units (EU)
Counts in Leaf Samples (LS)
Model expectations
Whiteflies
Predators
Whitefly nymphs
LS
EU
N
Ny-EU
N
Ad
(
O
)
N
N
t
N
O
x
N
L
i
N
M
m
Dead
Alive
Emerged
s
N
Ny
(
O
)
N
Ny
(
E
)
N
Ad
(
E
)
C
A
956
318
1
0
0
0
43
154
172
369
326
334
385
C
A
830
575
3
0
0
0
11
235
126
372
361
337
334
C
A
937
418
0
0
0
0
30
156
223
409
379
371
377
C
B
280
60
0
0
0
0
17
45
34
96
79
87
58
C
B
480
35
2
0
0
0
36
129
72
237
201
215
99
C
B
530
28
0
0
0
0
7
115
61
183
176
166
109
Li
A
670
198
0
0
2
0
38
221
70
329
291
298
205
Li
A
707
251
0
0
1
0
17
127
109
253
236
229
217
Li
A
911
289
0
0
1
0
17
132
140
289
272
262
279
Li
B
910
126
2
0
3
0
22
179
50
251
229
228
142
Li
B
350
44
1
0
1
0
16
101
38
155
139
141
55
Li
B
550
76
1
0
2
0
7
114
43
164
157
149
86
Li + Mm
A
405
233
2
0
2
7
33
124
85
242
209
200
149
Li + Mm
A
501
377
0
0
1
7
8
104
195
307
299
266
193
Li + Mm
A
300
102
0
0
1
2
48
60
46
154
106
114
99
Li + Mm
B
810
102
0
0
2
4
16
202
67
285
269
257
165
Li + Mm
B
150
22
0
0
0
9
48
15
13
76
28
30
13
Li + Mm
B
340
38
0
0
3
5
43
69
31
143
100
107
58
Li + Ox
A
1030
278
0
1
2
0
23
195
211
429
406
389
315
Li + Ox
A
717
293
0
1
0
0
48
185
134
367
319
333
220
Li + Ox
A
510
100
0
1
1
0
23
134
83
240
217
218
156
Li + Ox
B
30
5
0
0
0
0
0
9
4
13
13
12
5
Li + Ox
B
190
48
1
0
3
0
35
79
43
157
122
142
30
Li + Ox
B
530
86
1
0
1
0
20
107
37
164
144
149
83
Mm
A
406
107
0
0
0
5
19
107
45
171
152
148
156
Mm
A
368
257
0
0
1
5
24
64
84
172
148
146
139
Mm
A
865
455
0
0
0
7
32
159
142
333
301
301
348
Mm
B
380
60
0
0
0
5
36
107
44
187
151
155
71
Mm
B
70
10
2
0
0
7
33
12
1
46
13
15
5
Mm
B
290
41
0
0
0
8
39
59
41
139
100
105
50
Mm + Nt
A
540
166
3
0
0
3
118
115
97
330
212
242
176
Mm + Nt
A
717
214
5
0
0
6
26
136
83
245
219
190
247
Mm + Nt
A
340
142
6
0
0
3
87
38
69
194
107
128
100
Mm + Nt
B
260
37
7
0
0
1
105
34
39
178
73
72
24
Mm + Nt
B
370
20
0
0
0
2
81
52
46
179
98
87
41
Mm + Nt
B
520
98
12
0
0
6
106
94
33
233
127
132
67
Mm + Ox
A
406
302
0
2
0
5
64
177
82
323
259
277
155
Mm + Ox
A
513
150
1
1
0
8
62
86
82
230
168
203
202
Mm + Ox
A
1776
750
0
0
0
5
45
164
192
401
356
363
715
Mm + Ox
B
1230
244
8
2
0
8
25
172
61
258
233
234
253
Mm + Ox
B
250
24
0
1
0
5
38
68
49
155
117
106
39
Mm + Ox
B
280
77
5
2
0
12
40
56
45
141
101
102
46
Nt
A
213
210
3
0
0
0
31
66
67
164
133
127
73
Nt
A
466
143
11
0
0
0
45
105
48
198
153
161
168
Nt
A
390
170
11
0
0
0
59
115
143
317
258
254
139
Nt
B
130
4
13
0
0
0
62
11
1
74
12
28
11
Nt
B
320
67
15
0
0
0
94
25
34
153
59
72
34
Nt
B
230
24
0
0
0
0
70
44
39
153
83
66
22
Nt + Li
A
510
280
7
0
0
0
40
85
110
235
195
171
165
Nt + Li
A
685
263
12
0
4
0
94
91
125
310
216
238
233
Nt + Li
A
832
158
4
0
1
0
72
92
64
228
156
181
293
Nt + Li
B
400
18
7
0
1
0
113
19
10
142
29
74
47
Nt + Li
B
170
16
2
0
2
0
68
44
3
115
47
47
16
Nt + Li
B
430
39
2
0
0
0
67
97
46
210
143
112
52
Nt + Ox
A
271
135
6
3
0
0
80
59
65
204
124
119
70
Nt + Ox
A
1208
440
4
1
0
0
65
158
104
327
262
266
436
Nt + Ox
A
823
224
9
1
0
0
124
155
115
394
270
294
273
Nt + Ox
B
670
96
8
1
0
0
77
83
29
189
112
101
81
Nt + Ox
B
600
82
16
1
0
0
110
83
36
229
119
116
69
Nt + Ox
B
310
69
5
1
0
0
101
34
30
165
64
60
26
Ox
A
364
119
0
1
0
0
15
158
54
227
212
206
111
Ox
A
970
270
1
1
0
0
22
190
118
330
308
299
297
Ox
A
467
297
0
1
0
0
19
93
146
258
239
234
143
Ox
B
1010
95
3
0
0
0
12
215
60
287
275
260
158
Ox
B
410
25
1
0
0
0
22
95
67
184
162
167
64
Ox
B
640
70
1
1
0
0
6
91
69
166
160
150
100
Averages
A
655
263
45
128
110
283
238
239
229
B
427
60
46
82
39
167
121
121
71
N
Ny-EU: number of whitefly nymphs counted in each experimental unit; NAd(O): number of whitefly adults counted in each experimental unit; NNt, NOx, NLi, NMm: number of predator counted in each experimental unit for the species N. tenuis, O. lineatus, P. phrygianus, and M. melanotoma respectively; Dead: number of dead whitefly nymphs counted in leaf samples; Alive: number of alive whitefly nymphs counted in leaf samples; Emerged: number of pupal cases from which whitefly adult emerged counted in leaf samples; s: sample size (Dead + Alive + Emerged); NNy(O): total number of whitefly nymphs survived in leaf samples (Alive + Emerged); NNy(E) and NAd(E): expected numbers of whitefly nymphs and adults surviving in leaf samples and whole experimental unit respectively (estimated from model J9). In bold the three series of counts confronted with expectations to obtain maximum likelihood estimation of model parameters.


### 3.2. What Do Models Capture?

Table1 shows model selection results and parameter estimates coming from the analysis of data relative to whitefly nymph survival (NNy(O); models S1–S17). The best model selected, S14 (R2 = 0.9755):
(11)p(Ny)(C,i,B,ai,M);pOx′=pLi′=0;pMm′=pMm+Pj′;pNt′<pNt+Pj′,
specifies that there is no predation on nymphs from O. lineatus and P. phrygianus (pOx′=pLi′=0); there is some unknown mortality in all replicates (C); per capita predation rates of N. tenuis and M. melanotoma vary with predator identity (i); block effects are present and do not change with predator identity (B); whitefly nymph starting density affects per capita predation rates differently with predator identity (ai); per capita predation rates differ between single- and multiple-predator treatments (M). 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 (pNt′<pNt+Pj′;pNt′=0.1171,pNt+Pj′=0.1722). By contrast, M. melanotoma does not seem to change per capita predation rate when in combination with any other predator (pMm′=pMm+Pj′=0.3338). However, model S13 (wi=0.4473;ΔQAICc=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 (pMm′=0.2974;pMm+Li′=0.4972).Figure1 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 R2, the model provides a good description of the variation present in the data.Confrontation of number of whitefly nymphs survived in each leaf sample (NNy(O)) with expectations (N(Ny(E))±SD) provided by model S13.
(a)
(b)The other models presented in Table1 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 (ΔQAICc=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 (ai), see the simplification from S4 to S5. To evaluate strength of evidence relative to block effects on predation rates (B), see the simplification from S14 to S16 (ΔQAICc=25.7016). Also, the very high ΔQAICc 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.Table2 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 ΔQAICc, the information relative to predation on whitefly adults appears rather uncertain, and four models have similar weight (J8–J11). The best model selected, J9 (R2=0.1037 and R2=0.6243, based on the kNB of the NAd(O) and Emerged data series, resp.):
(12)p(Ny)(S13)ψ(B)p(Ad)(M);pNt′=pMm′=0;pLi′=pOx′=pLi+Ox′,
specifies that predation on nymphs is as in model S13; emergence rate of whitefly adults changes with blocks, ψ(B), being less likely in Block B, where average day temperature was lower (ψA=0.4443; ψB=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, p(Ad)(M); there is no evidence for predation on whitefly adults from N. tenuis and M. melanotoma (pNt′=pMm′=0); there is some evidence about predation in the Ox, Li, and Li + Ox treatments with similar predation rates (pLi′=pOx′=pLi+Ox′).Figure2 confronts the number of whitefly adults counted in each experimental unit with expectations and standard deviations (SDs) provided by best model, J9. By showing very high standard deviations and approximative matches between observations and expectations, the figure underlines how the model captures very little information in the adult counts data series as also suggested by the very low estimated R2 based on the kNB of the NAd(O) data series (0.1037).Confrontation of number of whitefly adults counted in each experimental unit (NAd(O)) with expectations (NAd(E) ± SD) provided by model J9.
(a)
(b)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 (pLi′=pOx′,J8); (2) spiders prey on whitefly adults in single-predator treatments and in combination, but not when they are with the mirids (pLi′=pOx′=pLi+Ox′=pOx+Li′,J10); there is no predation at all on whitefly adults, p(Ad)(k) (J11). Whereas models J8, J9, and J10 state similar things about spider predation on whitefly adults, it is important to note that support obtained from model J11 and the very low R2 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 (ΔQAICc 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ΔQAICc 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 R2 of the best model based on the kNB relative to the Emerged data series, which is rather high, 0.6243, relative to model ψ(B) versus ψ(·).


### 3.3. 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, Figure3 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.Figure 3
Model predictions about variation of per capita predation rates of mirids with species identity, whitefly nymph density and interactions between predators. Per capita predation rates ofM. melanotoma alone (pMm′) and in combination with P. phrygianus (pMm+Li′) decrease with increasing whitefly nymph density. Per capita predation rates of N. tenuis alone (pNt′) and in combination with any other predator (pNt+Pj′) tend to decrease less sharply with whitefly density. The two predators appear complementary across prey density. Parameter estimates taken from model J9.


### 3.4. Predation between Predators

As described in the methods section, all models ranked in Table1 assume that number of predatory Miridae changed during the experiment according to Table 4, thus, corresponding to the structure “MIGP: 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-MIGP), model “S14-MBM: different background mortality present, no IGP” has a ΔQAICc=5.5942, whereas model “S14-MNULL: no different mortality and no IGP” has a ΔQAICc=5.5978. The rather high ΔQAICc 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 ΔQAICc(5.5942) results from excluding IGP (S14-MBM), 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.


## 3.1. What Do Raw Data Suggest?

Table5 shows how counts of the 3 data series analysed varied among experimental units and between blocks: NNy(O) = 238±81 (SD) and 121±69 for blocks A and B, respectively; Emerged = 110±35 and 39±17 for A and B; NAd(O) = 263±142 and 60±42 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 (NNy-EU = 655±323 and 427±262; 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 NAd(O)/NNy-EU = 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 Table3 (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.Table 3
Mean mortality of whitefly nymphs (Dead/s) for each treatment and control (n=6).
Nt
Mm
Ox
Li
Nt
0.4187
Mm
0.4016
0.2562
Ox
0.4009
0.2009
0.0673
Li
0.4159
0.2438
0.1042
0.0792
C
0.0978
Dead: number of dead whitefly nymphs counted in leaf samples; s: sample size (Dead + Alive + Emerged); Nt, Mm, Ox, Li: predator treatments (N. tenuis, M. melanotoma, O. lineatus, and P. phrygianus resp.), C: control.Table 4
Mean number of predators recovered in each treatment (±SD;n=6).
Nt
Mm
Ox
Li
Nt
8.83
±
5.95
3.50
±
2.07
1.33
±
0.82
1.33
±
1.51
Mm
5.50
±
4.04
6.17
±
1.33
1.33
±
0.82
1.50
±
1.05
Ox
8.00
±
4.34
7.17
±
2.79
0.67
±
0.52
1.17
±
1.17
Li
5.67
±
3.83
5.67
±
2.50
0.50
±
0.55
1.67
±
0.82
Mean
7.00
5.63
0.96
1.42
Number of predators recovered for each experimental unit are given in Table5 (NNt, NMm, NOx, NLi), 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.Table 5
Number of whitefly nymphs, adults and predators counted at the end of the experiment and model expectations.
Treatment
Block
Counts in Experimental Units (EU)
Counts in Leaf Samples (LS)
Model expectations
Whiteflies
Predators
Whitefly nymphs
LS
EU
N
Ny-EU
N
Ad
(
O
)
N
N
t
N
O
x
N
L
i
N
M
m
Dead
Alive
Emerged
s
N
Ny
(
O
)
N
Ny
(
E
)
N
Ad
(
E
)
C
A
956
318
1
0
0
0
43
154
172
369
326
334
385
C
A
830
575
3
0
0
0
11
235
126
372
361
337
334
C
A
937
418
0
0
0
0
30
156
223
409
379
371
377
C
B
280
60
0
0
0
0
17
45
34
96
79
87
58
C
B
480
35
2
0
0
0
36
129
72
237
201
215
99
C
B
530
28
0
0
0
0
7
115
61
183
176
166
109
Li
A
670
198
0
0
2
0
38
221
70
329
291
298
205
Li
A
707
251
0
0
1
0
17
127
109
253
236
229
217
Li
A
911
289
0
0
1
0
17
132
140
289
272
262
279
Li
B
910
126
2
0
3
0
22
179
50
251
229
228
142
Li
B
350
44
1
0
1
0
16
101
38
155
139
141
55
Li
B
550
76
1
0
2
0
7
114
43
164
157
149
86
Li + Mm
A
405
233
2
0
2
7
33
124
85
242
209
200
149
Li + Mm
A
501
377
0
0
1
7
8
104
195
307
299
266
193
Li + Mm
A
300
102
0
0
1
2
48
60
46
154
106
114
99
Li + Mm
B
810
102
0
0
2
4
16
202
67
285
269
257
165
Li + Mm
B
150
22
0
0
0
9
48
15
13
76
28
30
13
Li + Mm
B
340
38
0
0
3
5
43
69
31
143
100
107
58
Li + Ox
A
1030
278
0
1
2
0
23
195
211
429
406
389
315
Li + Ox
A
717
293
0
1
0
0
48
185
134
367
319
333
220
Li + Ox
A
510
100
0
1
1
0
23
134
83
240
217
218
156
Li + Ox
B
30
5
0
0
0
0
0
9
4
13
13
12
5
Li + Ox
B
190
48
1
0
3
0
35
79
43
157
122
142
30
Li + Ox
B
530
86
1
0
1
0
20
107
37
164
144
149
83
Mm
A
406
107
0
0
0
5
19
107
45
171
152
148
156
Mm
A
368
257
0
0
1
5
24
64
84
172
148
146
139
Mm
A
865
455
0
0
0
7
32
159
142
333
301
301
348
Mm
B
380
60
0
0
0
5
36
107
44
187
151
155
71
Mm
B
70
10
2
0
0
7
33
12
1
46
13
15
5
Mm
B
290
41
0
0
0
8
39
59
41
139
100
105
50
Mm + Nt
A
540
166
3
0
0
3
118
115
97
330
212
242
176
Mm + Nt
A
717
214
5
0
0
6
26
136
83
245
219
190
247
Mm + Nt
A
340
142
6
0
0
3
87
38
69
194
107
128
100
Mm + Nt
B
260
37
7
0
0
1
105
34
39
178
73
72
24
Mm + Nt
B
370
20
0
0
0
2
81
52
46
179
98
87
41
Mm + Nt
B
520
98
12
0
0
6
106
94
33
233
127
132
67
Mm + Ox
A
406
302
0
2
0
5
64
177
82
323
259
277
155
Mm + Ox
A
513
150
1
1
0
8
62
86
82
230
168
203
202
Mm + Ox
A
1776
750
0
0
0
5
45
164
192
401
356
363
715
Mm + Ox
B
1230
244
8
2
0
8
25
172
61
258
233
234
253
Mm + Ox
B
250
24
0
1
0
5
38
68
49
155
117
106
39
Mm + Ox
B
280
77
5
2
0
12
40
56
45
141
101
102
46
Nt
A
213
210
3
0
0
0
31
66
67
164
133
127
73
Nt
A
466
143
11
0
0
0
45
105
48
198
153
161
168
Nt
A
390
170
11
0
0
0
59
115
143
317
258
254
139
Nt
B
130
4
13
0
0
0
62
11
1
74
12
28
11
Nt
B
320
67
15
0
0
0
94
25
34
153
59
72
34
Nt
B
230
24
0
0
0
0
70
44
39
153
83
66
22
Nt + Li
A
510
280
7
0
0
0
40
85
110
235
195
171
165
Nt + Li
A
685
263
12
0
4
0
94
91
125
310
216
238
233
Nt + Li
A
832
158
4
0
1
0
72
92
64
228
156
181
293
Nt + Li
B
400
18
7
0
1
0
113
19
10
142
29
74
47
Nt + Li
B
170
16
2
0
2
0
68
44
3
115
47
47
16
Nt + Li
B
430
39
2
0
0
0
67
97
46
210
143
112
52
Nt + Ox
A
271
135
6
3
0
0
80
59
65
204
124
119
70
Nt + Ox
A
1208
440
4
1
0
0
65
158
104
327
262
266
436
Nt + Ox
A
823
224
9
1
0
0
124
155
115
394
270
294
273
Nt + Ox
B
670
96
8
1
0
0
77
83
29
189
112
101
81
Nt + Ox
B
600
82
16
1
0
0
110
83
36
229
119
116
69
Nt + Ox
B
310
69
5
1
0
0
101
34
30
165
64
60
26
Ox
A
364
119
0
1
0
0
15
158
54
227
212
206
111
Ox
A
970
270
1
1
0
0
22
190
118
330
308
299
297
Ox
A
467
297
0
1
0
0
19
93
146
258
239
234
143
Ox
B
1010
95
3
0
0
0
12
215
60
287
275
260
158
Ox
B
410
25
1
0
0
0
22
95
67
184
162
167
64
Ox
B
640
70
1
1
0
0
6
91
69
166
160
150
100
Averages
A
655
263
45
128
110
283
238
239
229
B
427
60
46
82
39
167
121
121
71
N
Ny-EU: number of whitefly nymphs counted in each experimental unit; NAd(O): number of whitefly adults counted in each experimental unit; NNt, NOx, NLi, NMm: number of predator counted in each experimental unit for the species N. tenuis, O. lineatus, P. phrygianus, and M. melanotoma respectively; Dead: number of dead whitefly nymphs counted in leaf samples; Alive: number of alive whitefly nymphs counted in leaf samples; Emerged: number of pupal cases from which whitefly adult emerged counted in leaf samples; s: sample size (Dead + Alive + Emerged); NNy(O): total number of whitefly nymphs survived in leaf samples (Alive + Emerged); NNy(E) and NAd(E): expected numbers of whitefly nymphs and adults surviving in leaf samples and whole experimental unit respectively (estimated from model J9). In bold the three series of counts confronted with expectations to obtain maximum likelihood estimation of model parameters.


## 3.2. What Do Models Capture?

Table1 shows model selection results and parameter estimates coming from the analysis of data relative to whitefly nymph survival (NNy(O); models S1–S17). The best model selected, S14 (R2 = 0.9755):
(11)p(Ny)(C,i,B,ai,M);pOx′=pLi′=0;pMm′=pMm+Pj′;pNt′<pNt+Pj′,
specifies that there is no predation on nymphs from O. lineatus and P. phrygianus (pOx′=pLi′=0); there is some unknown mortality in all replicates (C); per capita predation rates of N. tenuis and M. melanotoma vary with predator identity (i); block effects are present and do not change with predator identity (B); whitefly nymph starting density affects per capita predation rates differently with predator identity (ai); per capita predation rates differ between single- and multiple-predator treatments (M). 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 (pNt′<pNt+Pj′;pNt′=0.1171,pNt+Pj′=0.1722). By contrast, M. melanotoma does not seem to change per capita predation rate when in combination with any other predator (pMm′=pMm+Pj′=0.3338). However, model S13 (wi=0.4473;ΔQAICc=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 (pMm′=0.2974;pMm+Li′=0.4972).Figure1 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 R2, the model provides a good description of the variation present in the data.Confrontation of number of whitefly nymphs survived in each leaf sample (NNy(O)) with expectations (N(Ny(E))±SD) provided by model S13.
(a)
(b)The other models presented in Table1 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 (ΔQAICc=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 (ai), see the simplification from S4 to S5. To evaluate strength of evidence relative to block effects on predation rates (B), see the simplification from S14 to S16 (ΔQAICc=25.7016). Also, the very high ΔQAICc 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.Table2 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 ΔQAICc, the information relative to predation on whitefly adults appears rather uncertain, and four models have similar weight (J8–J11). The best model selected, J9 (R2=0.1037 and R2=0.6243, based on the kNB of the NAd(O) and Emerged data series, resp.):
(12)p(Ny)(S13)ψ(B)p(Ad)(M);pNt′=pMm′=0;pLi′=pOx′=pLi+Ox′,
specifies that predation on nymphs is as in model S13; emergence rate of whitefly adults changes with blocks, ψ(B), being less likely in Block B, where average day temperature was lower (ψA=0.4443; ψB=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, p(Ad)(M); there is no evidence for predation on whitefly adults from N. tenuis and M. melanotoma (pNt′=pMm′=0); there is some evidence about predation in the Ox, Li, and Li + Ox treatments with similar predation rates (pLi′=pOx′=pLi+Ox′).Figure2 confronts the number of whitefly adults counted in each experimental unit with expectations and standard deviations (SDs) provided by best model, J9. By showing very high standard deviations and approximative matches between observations and expectations, the figure underlines how the model captures very little information in the adult counts data series as also suggested by the very low estimated R2 based on the kNB of the NAd(O) data series (0.1037).Confrontation of number of whitefly adults counted in each experimental unit (NAd(O)) with expectations (NAd(E) ± SD) provided by model J9.
(a)
(b)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 (pLi′=pOx′,J8); (2) spiders prey on whitefly adults in single-predator treatments and in combination, but not when they are with the mirids (pLi′=pOx′=pLi+Ox′=pOx+Li′,J10); there is no predation at all on whitefly adults, p(Ad)(k) (J11). Whereas models J8, J9, and J10 state similar things about spider predation on whitefly adults, it is important to note that support obtained from model J11 and the very low R2 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 (ΔQAICc 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ΔQAICc 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 R2 of the best model based on the kNB relative to the Emerged data series, which is rather high, 0.6243, relative to model ψ(B) versus ψ(·).


## 3.3. 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, Figure3 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.Figure 3
Model predictions about variation of per capita predation rates of mirids with species identity, whitefly nymph density and interactions between predators. Per capita predation rates ofM. melanotoma alone (pMm′) and in combination with P. phrygianus (pMm+Li′) decrease with increasing whitefly nymph density. Per capita predation rates of N. tenuis alone (pNt′) and in combination with any other predator (pNt+Pj′) tend to decrease less sharply with whitefly density. The two predators appear complementary across prey density. Parameter estimates taken from model J9.


## 3.4. Predation between Predators

As described in the methods section, all models ranked in Table1 assume that number of predatory Miridae changed during the experiment according to Table 4, thus, corresponding to the structure “MIGP: 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-MIGP), model “S14-MBM: different background mortality present, no IGP” has a ΔQAICc=5.5942, whereas model “S14-MNULL: no different mortality and no IGP” has a ΔQAICc=5.5978. The rather high ΔQAICc 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 ΔQAICc(5.5942) results from excluding IGP (S14-MBM), 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.


## 4. 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.


### 4.1. Different Predator Effectiveness

In general, predator-specific 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.


### 4.2. Contrasting Interactions between Predators

Ecological interactions between predators [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.On the opposite, prey risk enhancement is well supported by the constant increase of the estimated per capita predation rate ofN. tenuis when combined with any of the other predators (R2 and ΔQAICc strongly support this hypothesis; note also that this result holds with model S14-MNULL: 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.


### 4.3. 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. melanotomaappears 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 control with different prey densities. This complementarity should be considered in biological control strategies, where release of only one species of these two apparently redundant predators is generally suggested for whitefly control [61, 62].


### 4.4. Environmental Conditions Affect Stage Transition and Predation

Results showed block effects on emergence rates of whiteflies and on per capita predation rates of predatory Miridae on whitefly nymphs, likely due to the observed variation of environmental temperature between temporal blocks (Block A:27.7±0.78°C; Block B: 22.7±1.07°C). Emergence rates of whitefly adults from nymphs were reduced by the lower environmental temperature registered during Block B. By contrast, as the estimated additive factor βi that applies only to Block A is negative (see model 4 and 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 ai. 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).


### 4.5. 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 single-prey 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].


## 4.1. Different Predator Effectiveness

In general, predator-specific 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.


## 4.2. Contrasting Interactions between Predators

Ecological interactions between predators [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.On the opposite, prey risk enhancement is well supported by the constant increase of the estimated per capita predation rate ofN. tenuis when combined with any of the other predators (R2 and ΔQAICc strongly support this hypothesis; note also that this result holds with model S14-MNULL: 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.


## 4.3. 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. melanotomaappears 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 control with different prey densities. This complementarity should be considered in biological control strategies, where release of only one species of these two apparently redundant predators is generally suggested for whitefly control [61, 62].


## 4.4. Environmental Conditions Affect Stage Transition and Predation

Results showed block effects on emergence rates of whiteflies and on per capita predation rates of predatory Miridae on whitefly nymphs, likely due to the observed variation of environmental temperature between temporal blocks (Block A:27.7±0.78°C; Block B: 22.7±1.07°C). Emergence rates of whitefly adults from nymphs were reduced by the lower environmental temperature registered during Block B. By contrast, as the estimated additive factor βi that applies only to Block A is negative (see model 4 and 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 ai. 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).


## 4.5. 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 single-prey 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].


## 5. 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].



---

*Source: 619393-2013-10-20.xml*