Invertebrate community structure predicts natural pest control resilience to insecticide exposure

1. Biological pest control has become one of the central principles of ecological intensification in agriculture. However, invertebrate natural enemies within agri cultural ecosystems are exposed to a myriad of different pesticides at both lethal and sub-lethal doses, that may limit their capacity to carry out pest control. An important question is how underlying diversity in invertebrate predator species, linked to their unique susceptibility to insecticides, can act to increase the resil ience of natural pest control. 2. We explore this issue by assessing the effects of sub-lethal insecticide exposure on the predation rates of 12 generalist predators feeding on the aphid Sitobion avenae (Aphididae). Predation rates within a 24-hr period were assessed (preda

important question is how underlying diversity in invertebrate predator species, linked to their unique susceptibility to insecticides, can act to increase the resilience of natural pest control.
2. We explore this issue by assessing the effects of sub-lethal insecticide exposure on the predation rates of 12 generalist predators feeding on the aphid Sitobion avenae (Aphididae). Predation rates within a 24-hr period were assessed (predation assessment) for each species after receiving one of the following treatments: (a) no prior deltamethrin exposure before the predation assessment (control); (b) deltamethrin exposure immediately before the predation assessment (resistance) and (c) deltamethrin exposure 5 days before the predation assessment (recovery).
Extrapolating from these species-specific measures of resistance and recovery, we predicted the resilience of community level predation to insecticide exposure for predator communities associated with 256 arable fields in the UK.
3. There was large variation in sub-lethal effects of the insecticide between even closely related species. This ranged from species showing no change in predation rates following sub-lethal insecticide exposure (high resistance), species showing only immediate depressed feeding rates after 24 hr (high recovery) or those with depressed feeding rates after 5 days (low resistance and recovery).
4. The community level analysis showed that resistance and recovery of natural pest control was predicted by both community phylogenetic diversity (positively) and weighted mean body mass (negatively). However, the removal of numerically dominant species from the analysis modified these effects.

Synthesis and applications.
Our results highlight the role of community diversity in maintaining the resilience of natural pest control following insecticide use.
Importantly, less diverse assemblages dominated by predator species that show

| INTRODUC TI ON
The utilisation of biological pest control in agricultural ecosystems has become central to the concept of ecological intensification, whereby farming systems integrate natural ecosystem services to offset anthropogenic inputs (Bommarco, Kleijn, & Potts, 2013).
There is strong evidence to suggest that natural predation can be optimised in combination with conventional agro-chemical control methods, with the potential to support integrated pest management strategies (Naranjo & Ellsworth, 2009). For example, Naranjo & Ellsworth, 2009 showed that multiple applications of a broad-spectrum insecticide (which strongly depressed natural enemy populations) were needed to control Bemisia tabaci in cotton production, compared with a single application of insect growth regulator, which due to its mode of action, had less impact on natural predator populations. This approach maximised natural pest control, providing the same level of control as broad-spectrum insecticides and saved growers up to $200 million (Naranjo & Ellsworth, 2009).
The effectiveness of integrated pest management will be maximised where one part of the control strategy (e.g. insecticides) does not degrade the other (e.g. beneficial predators). In many instances, broad-spectrum insecticides not only diminish withinfield natural enemy populations, but nullify attempts to increase their populations (Gagic et al., 2019;Ricci et al., 2019). The obvious mechanism for the negative effects of insecticides on natural enemies is exposure leading to direct mortality (Guedes, Smagghe, Stark, & Desneux, 2016;Stark & Banks, 2003). Historically, ecotoxicological testing has focused on median lethal dose (LD 50 ) or lethal concentration (LC 50 ) values necessary to kill 50% of a population (Stark & Banks, 2001. However, direct mortality is only one outcome of exposure to insecticides. While useful for describing the immediate toxicity of a chemical, insecticides persist in the environment for varying time periods after application (Goulson, 2013;Tang et al., 2018). This can result in longer-term exposure at sub-lethal doses that can impact on the biological viability of populations via other mechanisms, such as low fecundity (Desneux, Decourtye, & Delpuech, 2007). Such sub-lethal doses can also affect behaviours that could impact on ecosystem service delivery. For example, sub-lethal doses of pyrethroid and organophosphorus insecticides can impair locomotion of spiders and beetles for up to 9 days (Baatrup & Bayley, 1993;Singh, Port, & Walters, 2001;Tooming, Merivee, Must, Sibul, & Williams, 2014).
Many studies of sub-lethal effects are at the level of the individual, which is valuable for determining the range responses for an insecticide, however it is difficult to extrapolate these and apply them at the community level.
Understanding sub-lethal effects of insecticides on predators at the community scale will help to determine how resilient pest control ecosystem services are under pest management strategies that include chemical control. Resilience is fundamental to providing stable ecosystem service delivery, and can be broken down into two components; the first is resistance, which in the context of natural pest control describes how much predation deviates compared with baseline levels following exposure. The second is recovery, which can be viewed as the ability of pest control to return baseline levels a time after exposure (Beller et al., 2019;Kohler et al., 2017). The interplay between resistance and recovery within communities of natural enemies following exposure to insecticides will determine the efficacy of integrated pest management strategies. Additionally, as natural pest control is underpinned by components of community structure, including functional diversity, the ability of biodiversity to increase resistance and recovery is of considerable importance (Greenop, Woodcock, Wilby, Cook, & Pywell, 2018).
A key challenge remains in bridging the gap between responses of individual predators to insecticides in the laboratory and how this impacts the resilience of pest control services in real-world agricultural systems. Additionally, it is also important to understand if components of community diversity could help mitigate negative effects of insecticides on ecosystem services. In this study, we combine a laboratory experiment with data from 256 invertebrate community samples from UK arable fields collected as part of the farm scale evaluation (FSE; Firbank et al., 2003). We assess the predation responses of 12 common predators of the grain aphid Sitobion avenae under different exposures to the historically widely used insecticide deltamethrin. We model these effects based on real-world predator communities to identify how different components of diversity mediate effects of insecticide exposure on predation. We focus on generalist predators due to their importance as biocontrol agents in agricultural ecosystems (Symondson, Sunderland, & Greenstone, 2002). We address the following predictions, (a) predators will show a decrease in predation in the 24 hr immediately following sub-lethal exposure to an insecticide, but will demonstrate partial recovery after 5 days (Baatrup & Bayley, 1993;Tooming et al., 2014); (b) at the community level, greater diversity will increase the resistance and recovery of predation in response to insecticide exposure, due to mechanisms including the insurance effect whereby there is an increased likelihood that a resilient predator will be present in more diverse assemblages (Oliver et al., 2015).

| Study species
We determined the effects of a field-realistic exposure of the insecticide deltamethrin on 11 species of generalist predator and two species of ladybird, which are predominantly specialist aphid predators, with the aim of quantifying both the resistance and recovery in their feeding rates on aphids. We assessed predation of 10 species of ground beetles ( Table S1). Abax parallelepipedus was not found to feed on the aphids in any of the experiments. The species were caught in pitfall traps or hand collected in Oxfordshire between May and July 2019 and were kept in a controlled temperature room at 18°C (16 hr L:8 hr D cycle). Predators were kept individually in Petri dishes including moist tissue for a maximum of 7 days before the start of the experiment and fed with flightless drosophila, rehydrated mealworm and aphids.

| Assessing resistance and recovery in predation rates
For all beetle species, we assessed predation on the grain aphid S. avenae, an important aphid pest of cereals frequently used as a model species for measuring pest control services (Bosem Baillod, Tscharntke, Clough, & Batáry, 2017;Mansion-Vaquié, Ferrante, Cook, Pell, & Lövei, 2017). We wanted to determine the ability of predators to predate on the pest species S. avenae within a 24-hr period following insecticide exposure. Deltamethrin was chosen as a historically widely used broad-spectrum insecticide representative of the pyrethroid class (applied to 54,112 ha of arable cropland in the UK in 2018; Garthwaite et al., 2018). We do not propose that the responses to this insecticide will be representative of all insecticides, rather this provides a baseline for understanding the breadth of species responses. Based on field estimations of ladybird exposure rates undertaken by Wiles and Jepson, (1994) 3.1 ng a.i indiv. −1 was used as a standard exposure rate within the study.
Prior to the predation assessments, predators were starved for 5 days during which time they were exposed to one of the three following insecticide treatments: (a) no prior deltamethrin exposure before the predation assessment (control); (b) deltamethrin exposure immediately before the predation assessment i.e. exposure at the end of the starvation period (resistance) and (c) deltamethrin exposure 5 days before the predation assessment i.e. exposure at the beginning of the starvation period (recovery). Following Everts et al. (1991), each individual was treated with either water or the deltamethrin treatment. The water only control was used to account for the potential effects of liquid application independent of the effects deltamethrin might have on the predators. Following this protocol, at the start of the starvation period, individuals in the recovery treatment received 3.1 ng of deltamethrin dissolved in 1 μl water, while the other treatments received 1 μl of water applied to the dorsal side of the abdomen using a micropipette. Then after the 5-day starvation period individuals in the resistance treatment received 3.1 ng of deltamethrin dissolved in 1 μl water, while predators in the control and recovery treatments received 1 μl of water.
This approach meant all treatments could be carried out at the same time within a block. In all cases the application of deltamethrin occurred at approximately 12:00 hr. After the starvation period predators were weighed and introduced into opaque plastic arenas (L = 220 mm × W = 155 mm × H = 150 mm) with sides that were coated in Fluon ® (AGC) a synthetic fluropolymer that was used to stop aphids climbing up the side of the arena (Hentley et al., 2016).
Each arena contained 20 adult S. avenae aphids on a piece of wheat leaf 2-cm long and was lined with moist paper towel to provide moisture and habitat. Predators were given 24 hr to feed on the aphids, after which the predator was removed and weighed, and all the adult aphids were then counted. Predation was only assessed if the individual was alive at the end of the predation assessment. All deaths were recorded from the beginning of the starvation period until the end of the predation assessment for each block. This was done to determine whether deltamethrin exposure increased mortality during the starvation period.
Experiments were carried out in 28 blocks from May to July 2019 in a controlled environment room kept at 18°C (16 hr L:8 hr D cycle).
Species were tested based on their availability within blocks, where possible at least one replicate for each treatment for each species was carried out at the same time. We include a random effect to account for differences within species between blocks (see Section 2.2.1). For each predator species we obtained between eight and 10 replicates for each treatment using a new individual for every replicate (total replicates for each species are given in Appendix S1: Table S2). For A. dorsalis and B. bullatus we were only able to catch enough individuals to carry out the control and resistance exposure treatment. We also carried out 10 control replicates without predators to determine if there was a loss of aphids for reasons other than predation. Within the 24-hr assessment period there was no loss of aphids in this control.

| Statistics (part 1)
To determine the effects of the deltamethrin insecticide treatments on predation rates, we fitted Bayesian generalised linear mixed models to each predator species using the brms package in RStudio (Bürkner, 2017;Carpenter et al., 2017;R Core Team, 2020). The response variable was the proportion of aphids eaten and the explanatory variable Insecticide treatment (three levels: Control, Resistance and Recovery). All models included a temporal block descriptor as a random effect. Depending on the responses of individual species, models were fitted either using: (a) a binomial model; (b) a binomial model with an observation level random effect to account for overdispersion; or (c) a beta-binomial model to account for overdispersion, all with a logit link function. Model selection was based on which approach best addressed overdispersion using either k-fold (10 folds) or leave-one-out (loo) validation (Harrison, 2015;Vehtari, Gelman, & Gabry, 2017). Following prior sensitivity analysis (see Appendix S2), we used a Normal (μ = 0, σ = 2.5) prior on all main effects. Models were run with four chains for 4,000 iterations with 1,000 burn in iterations. Fit was based on posterior predictive checks, Rhat values <1.05 and inspection of residual plots (Gelman & Rubin, 1992). We calculated the mean posterior distribution of differences between treatment levels and 95% credible interval (CI) to determine the effects of deltamethrin on each predator species. Where credible intervals do not include zero indicates a significant effect. All results are given on the log odds scale. To provide context a log odds of 0 is equal to a probability of 0.5 (i.e. 50% of aphids consumed).

| Community resistance and recovery
To determine the extent to which insecticide exposure could impact pest control at the community level, we extrapolated responses found for individual species in the laboratory in terms of their resistance and recovery. This was undertaken for 256 real arable farm communities recorded as part of the UK FSE (Firbank et al., 2003). The FSE was setup to test the impacts of herbicide tolerant crops on agricultural wildlife and included wide scale monitoring of biodiversity across the following crop types: sugar beet (25.78% of fields), maize (22.66%), spring oilseed rape (26.17%) and winter oilseed rape (25.39%; Firbank et al., 2003). We selected 10 predator species from those tested in the laboratory experiment to calculate community responses at the FSE sites representing a mean of 47.81% (SD = 21.97%) of all predator species abundances present across the FSE data (see Appendix S3 for full details). Only 10 species were used as the predator H. axyridis was not present at any of the sites, possibly because the FSE trials took place before this invasive species first appeared in the UK (Majerus, Strawson, & Roy, 2006). Philonthus cognatus was excluded from FSE analysis as staphlynids were not identified to species level in the FSE data and A. parallelepipedus was excluded as it did not predate on aphids in the first experiment. Note, that for the recovery treatment we did not have data for A. dorsalis and B. bullatus due to a lack of captured individuals. For each site we calculated two response ratios giving a measure of community resistance and recovery based on the abundances of the 10 species (i.e. those we studied in the feeding experiment) at the FSE sites and their predation responses derived from the feeding experiment, these were: (a) the ratio of change in predation within a 24-hr period immediately after insecticide exposure (resistance; control estimated predation at a site/resistance estimated predation at a site) and (b) the ratio of change in predation 5 days after exposure (recovery; control estimated predation at a site/recovery estimated predation at a site; Appendix S3 for full methodological details). As random sampling was involved in the calculation of the response ratios they were derived 100 times for each site. For the resistance metric this gave 100 datasets each consisting of response ratios for 255 sites (n = 255) and for the recovery metric this gave response ratios for 254 sites (n = 254; sites were removed where they contained only a single species we had feeding data for, which meant phylogenetic diversity could not be calculated). While we focus here on sub-lethal effects, a certain number of individuals did die across all treatments during the experiment. To account for any potential differences between treatments in mortality, we carried out the modelling first including only sub-lethal effects, and then factoring in mortality by multiplying the abundance at each site for each species by its probability of survival derived from the laboratory experiment. This was also repeated 100 times for each site.
The response ratios were then used to determine how components of community diversity could mitigate effects of insecticide exposure on pest control. We used the following explanatory variables describing invertebrate community structure at each of the FSE sites: abundance (total number of individuals); species richness (count of number of species); and community evenness (Pielou's measure of species evenness Smith & Wilson, 1996), which was removed due to intercorrelations with functional diversity and community weighted mean (CWM) body mass (see Appendix S4 for correlation matrices). All measures have been linked to community resilience (Feit, Blüthgen, Traugott, & Jonsson, 2019;Loreau & de Mazancourt, 2013;Oliver et al., 2015). We also considered four metrics describing the functional trait structure of the communities, these were: (a) CWM body mass, where body mass has been shown to mitigate toxicity to insecticides (Wiles & Jepson, 1992); (b) CWM flight capacity (macropterous, brachypterous or dimorphic), as the process of opening the chitinous wing cases was considered to be a factor increasing exposure risk. None of the wing type CWMs were included in the final models due to intercorrelations with other variables (Appendix S4); (c) Functional diversity, represented by functional dispersion (an abundance weighted measure of functional diversity) was calculated using both body mass and wing type in the fd package (Laliberte & Legendre, 2010) and (d) Phylogenetic diversity, used to predict the potential for intrinsic differences in sensitivity to insecticides based on phylogenetic history related to toxicokinetic and toxicodynamic processes (Rubach et al., 2011). Phylogenetic diversity was abundance-weighted and was derived using the mean pairwise taxonomic relatedness of a taxonomy surrogate, standardised to correct for intercorrelations with species richness (Order, Family, Sub-family, Tribe, Genus and species). Phylogenetic diversity was assessed in the picante package in r using ses.mpd function (Kembel et al., 2010). Traits and taxonomy of species are included in Appendix S1: Table S2.

| Statistics (part 2)
We fitted a Bayesian linear model to each of the 100 generated datasets using rstanarm package in r (Goodrich, Gabry, Ali, & Brilleman, 2018). Each model was fitted with the explanatory variables described above (reference model). All models used a gaussian response distribution with either the resistance or recovery log response ratio as the response variable. Models were fitted using weakly informative Normal (0, 10) prior on the intercept, and a regularised horseshoe prior on the fixed effects (Goodrich et al., 2018;Piironen & Vehtari, 2017a). All models were run with four chains for 3,000 iterations and 1,000 warmup iterations (Goodrich et al., 2018). Fit was based on posterior predictive checks, Rhat values <1.05 and inspection of residual plots (Gelman & Rubin, 1992). From this starting point we then carried out projective predictive model selection on the reference model to determine a subset of parameters that best predicted community predation responses to deltamethrin exposure without an increase in predictive error (Piironen, Paasiniemi, & Vehtari, 2018;Piironen & Vehtari, 2017b). Variable selection was carried out using the cv_varsel function validated by 10-fold cross validation in the projpred package (Piironen et al., 2018;Vehtari et al., 2017). We present the percentage inclusion of all predictor  Table S1).

| Sensitivity analysis
We also conducted a sensitivity analysis for the community measures of resistance and recovery due to the relatively small number of species tested in order to determine how sensitive our findings were to the absence of each species. To test sensitivity we sequentially removed each species and then recalculated all our diversity metrics and response ratios using the same modelling procedure described above. We present all results from the sensitivity analyses in Appendix S6.

| Resistance and recovery of individual predator species
Of the predators sampled, the ladybird H. axyridis had the highest predation in the control treatment and the lowest predation observed was by N. brevicollis (Figure 1a). The four most efficacious predator species in the control showed reductions in feeding in the resistance treatment. The recovery treatment showed that for these species, feeding rates returned to levels indistinguishable from the control (excluding A. dorsalis where recovery was not assessed), however, survival for the recovery treatment was lower for H. axyridis and P. cognatus compared to the control (Figures 1b and 2; Appendix S5: Table S2). In contrast to H. axyridis, C. septempunctata predation was unaffected by exposure to deltamethrin, and while it still suffered mortality in the pesticide exposure treatments, this was lower than that observed in H. axyridis (Figure 2; Appendix S5: Table S2). Similarly, for both the Harpalus species, A. plebja, B. bullatus and N. brevicollis predation was found not to be strongly affected by exposure to sub-lethal levels of deltamethrin (Figure 1b). Pterostichus madidus showed poor recovery with a depression in feeding rates in response to deltamethrin compared to the control that persisted for 5 days  Table S2).

Resistance
Focusing on sub-lethal effects, all the highest performing models describing the resistance of communities to deltamethrin included phylogenetic diversity and CWM body mass (Table 1). Functional diversity was included in 1% of models, and species richness and abundance in 6% of models (Table 1) (Table 1). Functional diversity, abundance and species richness were again only included in a subset of models with coefficients that had credible intervals that overlapped zero (Table 1)

Resistance sensitivity
The removal of three species P. melanarius, P. madidus and N. brevicollis strongly affected the impact of the community diversity measures on resistance. Notably, when either P. melanarius or P. madidus were removed from both the sub-lethal and mortality analyses phylogenetic diversity no longer had a clear effect on resistance (see Appendix S6: Table S1 and Table S2). The coefficients for phylogenetic diversity were negative, however the credible intervals overlapped zero (see Appendix S6: Table S1 and Table S2). Where N. brevicollis was removed from the sub-lethal analysis CWM body mass had a negative effect on resistance but the credible intervals overlapped zero, whereas abundance had a clear negative effect (Appendix S6: Table   S1). The removal of N. brevicollis from the mortality analysis caused CWM body mass, abundance and phylogenetic diversity all to have credible intervals that overlapped zero (Appendix S6; Table S2).

Recovery
The results suggested that while recovery had occurred after 5 days there was still evidence of a depression in natural pest con-  (Table 1). Species richness and abundance were included in 99% and 58% of models, respectively, and functional diversity in 32%. Communities with greater phylogenetic diversity and species richness showed higher recovery (Table 1; Figure 4a,d). In comparison, CWM body mass and abundance decreased recovery to deltamethrin (Table 1; F I G U R E 1 (a) The control (no exposure to deltamethrin) log odds ratio of aphid predation for 12 predator species in a 24-hr period (predation assessment) analysed using Bayesian mixed models; (b) resistance and recovery show the difference in the log odds of aphid predation compared with the control, where individuals were exposed to a sublethal dose of deltamethrin immediately before the predation assessment (resistance; blue) and 5 days before the predation assessment (recovery; blue). Points are means and error bars show lower and upper 95% credible intervals. The dashed line indicates no difference from the control treatment F I G U R E 2 Overall survival of each species across the experimental period. Treatments were: no exposure to deltamethrin prior to a 24-hr predation assessment on the aphid species Sitobion avenae (control); exposure to a sub-lethal dose of deltamethrin immediately before the predation assessment (resistance); and exposure to the same dose of deltamethrin 5 days before the predation assessment (recovery) TA B L E 1 The minimum and maximum coefficient and percentage inclusion for each variable included in the highest performing Bayesian sub-models. Resistance refers to the log response ratio that estimated change in predation (compared with unexposed communities) following a sub-lethal dose of deltamethrin within a 24-hr period. Recovery refers to the log response ratio that estimated the change in predation 5 days after exposure to the same dose. Bold parameters and values indicate credible intervals for both the minimum and maximum estimate do not overlap zero  Figure 4d). As was seen in the assessment for resistance, including mortality effects did not qualitatively alter model predictions. Both phylogenetic diversity and species richness had a positive effect and were included in 100% of models (Table 1).
CWM body mass had a larger effect compared with models only including sub-lethal effects and was included in 100% of models (Table 1). Abundance was included in 29% of models and had a negative effect on the response ratio, as did functional diversity that was included in 10% of models, however had 95% credible intervals that overlapped zero (Table 1).

Recovery sensitivity
Similar to the resistance sensitivity analyses, the removal P. melanarius and P. madidus had a strong effect on the effects of the diversity

| Individual species resistance and recovery
The most common approach to assessing the toxicity of pesticides is to use representative species and a measurement of the dose required to kill 50% of the population (LD 50 or LC 50 ; Desneux et al., 2007). However, our results support other research (Desneux et al., 2007;Everts et al., 1991;Wiles & Jepson, 1994) that suggest sub-lethal exposure to insecticides, at doses below LD 50 or LC 50 values, can have impacts on the predation capacity of generalist predators which could affect natural pest control ecosystem services. The arthropod predators investigated were characterised by large variation in both resistance and recovery following sublethal exposure to deltamethrin, providing mixed support for our first prediction. This variation in terms of individual species resistance and recovery is a product of their unique toxicokinetics and toxicodynamics (Rubach et al., 2011). The utilisation of biomarker approaches may offer the ability in the future to identify the mechanistic differences that occur between different taxa, not in terms of just lethal but also sub-lethal effects (Desneux et al., 2007).
However, even the use of biomarker approaches has demonstrated large differences between closely related species in their mechanisms for dealing with toxicants (Spurgeon, Svendsen, Rimmer, Hopkin, & Weeks, 2000;Trekels, Van de Meutter, Bervoets, & Stoks, 2012). Addressing why species susceptibility to pesticides shows high variability will prove an important step in predicting how novel pesticides could impact communities in agricultural fields (Guedes et al., 2016).

| Community level resistance and recovery
When all species were included in the analyses, we found evidence that greater phylogenetic diversity and species richness to a certain extent, can increase the resistance and recovery of pest control ecosystem services, providing support for the insurance hypothesis (Balvanera et al., 2006). Theoretically, individuals in phylogenetically diverse communities may be less likely to share similar mechanisms for dealing with toxicants (see Guénard, von der Ohe, Walker, Lek, & Legendre, 2014).

F I G U R E 4
The impact of the diversity measures on the recovery log response ratio of change in predation following exposure to sub-lethal doses of deltamethrin, compared with unexposed communities. Solid line shows the mean and shaded areas the 95% credible intervals. All other variables included in the models were held at their mean.
(a) phylogenetic diversity, (b) community weighted mean body mass, (c) functional diversity, (d) species richness and (e) log abundance Therefore, where phylogenetic diversity is greater there is an increased chance of a species which shows high resistance or recovery to the insecticide and is therefore able to maintain predation (Balvanera et al., 2006;Oliver et al., 2015). However, the positive effects of phylogenetic diversity, and in certain cases species richness, were dependent on the inclusion of P. madidus and P. melanarius in the analysis. The species-dependent beneficial effects of phylogenetic diversity and species richness are indicative of a sampling effect where greater diversity increases the chance of a species that has a strong impact on ecosystem functioning being present within a community (Cadotte, Dinnage, & Tilman, 2012;Davies, Urban, Rayfield, Cadotte, & Peres-Neto, 2016;Gravel et al., 2012;Tilman, Lehman, & Thomson, 1997

| Caveats
This study does not consider the relative magnitude of predation, but rather focuses on the expected relative change in predation of communities when exposed to environmental stress. Our results do not lead to the conclusion that communities with greater phylogenetic diversity will provide a greater magnitude of predation than less diverse systems, rather that these systems are estimated to be more resistant and recover faster in their capacity to provide pest control when exposed to insecticides. We focus on a limited subset of predators with our results describing the response of a mean of 47.81% (SD = 21.97%) of the community depending on the field in the FSE data. Therefore, more research is needed to deter-

| CON CLUS IONS
Overall, we have found evidence to suggest that increasing phylogenetic diversity in agricultural ecosystems could increase the resistance and recovery of pest control ecosystem services following insecticide applications. Our results also highlight that these effects are likely to be strongly driven by the dominant predators within the community. In this study, we only consider a low dose of a single chemical, while in typical agricultural systems predators are exposed to a multitude of different pesticides. Therefore, the extent to which community diversity is likely to increase the resilience of pest control ecosystem services under these far more complicated conditions is difficult to predict. It seems reasonable to suggest that increasing the dose of broad-spectrum insecticides, such as deltamethrin, would further dampen both the resistance and recovery of natural pest control (Desneux, Rafalimanana, & Kaiser, 2004;Gyldenkaerne, Ravn, & Halling-Sørensen, 2000).
Current evidence from field studies would suggest that present management strategies (increasing habitat complexity or implementing field margins) may fail to strongly promote resilience of natural pest control, at least in response to commonly used insecticides (Gagic et al., 2019;Ricci et al., 2019). In order for integrated pest management strategies to become a viable option in open arable systems, determining chemicals that maximally impact the pest species while having minimal impact on beneficial invertebrates is an important step. Furthermore, integrating chemical control with land management that has the goal of increasing biodiversity, would appear to be fundamental in ensuring that ecological intensification strategies (such as the utilisation of natural pest control) are effective and deliver the desired outcomes (Gagic et al., 2019;Ricci et al., 2019).

ACK N OWLED G EM ENTS
The research was funded by the Natural Environment Research Council (NERC) under research programme NE/N018125/1 ASSIST -Achieving Sustainable Agricultural Systems www.assist.ceh.ac.uk.
ASSIST is an initiative jointly supported by NERC and the Biotechnology and Biological Sciences Research Council (BBSRC). B.A.W. is also funded under CHEMPOP (NEC06550) and CHEMMIX (NEC06567).

AUTH O R S ' CO NTR I B UTI O N S
All authors conceived the ideas and designed the methodology; A.G.
collected the data, analysed the data and led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data available from the Dryad Digital Repository https://doi. org/10.5061/dryad.66t1g 1k0f (Greenop, Cook, Wilby, Pywell, & Woodcock, 2020).