DOI: 10.5061/dryad.xksn02vhw
title: Bird predation and landscape context shape arthropod communities on broccoli

description: Birds increase crop yields via consumption of pests in some contexts but
 disrupt pest control via intraguild predation in others. Landscape
 complexity acts as an inconsistent mediator, sometimes increasing,
 decreasing, or not impacting pest control. Here, we examined how landscape
 context and seasonal variation mediate the impact of birds on arthropod
 pests and natural enemies, leaf damage, and yields of broccoli (Brassica
 oleracea) on highly diversified farms that spanned the USA West Coast. Our
 study had two complementary components: a bird exclusion experiment and
 molecular diet analysis of 357 fecal samples collected from the most
 commonly captured bird species that also foraged in Brassica fields
 (American Goldfinch, Spinus tristis; American Robin, Turdus migratorius;
 Savannah Sparrow, Passerculus sandwichensis; Song Sparrow, Melospiza
 melodia; and White-crowned Sparrow, Zonotrichia leucophrys). Bird access
 yielded higher, rather than lower, numbers of pest aphids and increased
 their parasitism, while no other arthropods examined were consistently
 impacted. Independent of bird presence, percent natural cover in the
 landscape sometimes increased and sometimes decreased densities of
 arthropods in the mid-growth period, with diminishing impacts in the
 late-growth period. Herbivore feeding damage to broccoli leaves decreased
 with increasing amounts of natural land cover and in the late-growth
 period. Molecular diet analysis revealed that Brassica pests and predatory
 arthropods were relatively uncommon prey for birds. Landscape context did
 not alter the prey items found in bird diets. Altogether, our
 bird-exclusion experiment and molecular diet analysis suggested that birds
 have relatively modest impacts on the arthropods associated with broccoli
 plantings. More broadly, the limited support in our study for net natural
 pest control services suggests that financial incentives may be required
 to encourage the adoption of bird friendly farming practices in certain
 cropping systems.
description: <strong>Study Sites</strong>
 Across three years (2016 – 2018), we conducted surveys to
 elucidate the net effects of birds on diversified crop production on farms
 throughout Washington, Oregon, and California, USA, previously described
 in Smith et al. (2020a,b, 2021). Farms grew a range of crops including
 cereals (e.g., corn, wheat, barley), vegetables and melons (e.g.,
 brassicas, leafy vegetables), fruits and nuts (e.g., citrus fruits,
 grapes, berries, walnuts), oilseed crops (e.g., olives, sunflower), roots
 (e.g., potatoes), spice crops (e.g., chilies, peppers, fennel), beverage
 crops (e.g., tea), medicinal crops, commercial flowers, and grasses and
 fodder crops, among others. All farms grew brassica crops including
 broccoli and/or kale. Farms in the study were generally small scale (mean
 = 25.3 ha ± 47.2 (SE); range = 0.44 – 272.2 ha). We conducted mist-netting
 to obtain avian fecal samples for molecular diet analysis in 2016 and 2017
 on 47 farms (Figure 1). In 2018, we recruited a subset consisting of 13
 cooperating farms from Washington and Oregon to participate in our bird
 exclusion experiment. We asked all farms in the Pacific Northwest region
 to participate, and of those, 13 agreed. We limited our experiment to the
 Pacific Northwest (i.e. did not conduct the experiment in California) for
 feasibility. <strong>Landscape
 context</strong> To characterize landscape
 context, we calculated the percent natural/semi-natural (hereafter,
 “natural”) cover based on the 2016 National Land Cover Database (Homer et
 al. 2012) using a 1000 m radius buffer centered at the cage location for
 our bird exclusion experiment and at the farm center for molecular diet
 analyses using Program R and FRAGSTATS 4.1 (McGarigal and Marks 1994).
 Natural cover included forest (deciduous, evergreen, and mixed), scrubland
 (dwarf scrub and shrub/scrub), herbaceous (grassland/herbaceous,
 sedge/herbaceous, lichens, and moss), and wetland categories (woody and
 emergent herbaceous wetlands). Categories not included in natural cover
 were water, ice/snow, developed, barren, pasture/hay, and cultivated crop
 classes. We used a 1000 m buffer to encompass the daily activity centers
 for both arthropods and birds and because it is commonly used for
 arthropod assessments in agricultural systems (Lichtenberg et al. 2017,
 Billerman et al. 2020, Tamburini et al. 2020).
  <strong>Bird Exclusion
 Experiments</strong> At each of the 13 farms in
 our exclusion experiment, we established 3 bird exclosures and 3 controls
 placed consecutively in a randomized order, each surrounding 4 broccoli
 plants (Figure 2; Martin et al. 2013, Maas et al. 2019). In all cases, our
 experimental broccoli plantings were established within broccoli fields
 that had been planted by the cooperating growers between 28 June and 23
 July 2018 (Figure 2D). We installed cages at the edges of plantings
 closest to seminatural elements to control for possible variation in pest
 control services due to edge effects (Olimpi et al. 2020). Cages were
 installed between the day of planting and 6 days post-planting. Cages were
 built from PVC pipe and covered in #6 monofilament 3.81-cm mesh netting
 (Figure 2C; Memphis Net & Twine, Memphis, TN, USA). Exclosures
 were 0.61-m tall and individualized in length and width to fit the crop
 spacing each grower used such that the cage would enclose 4 broccoli
 plants (range = 0.48 x 0.79 m to 1.2 x 1.2 m). Control plots were marked
 by thin twine strung on thin posts (~2-mm diameter) to discourage perching
 by birds (Maas et al. 2019). We conducted
 mid-broccoli-growth (9 Aug 18 to 17 Aug 18) and late-broccoli-growth (7
 Sep 18 to 23 Sep 18) surveys of arthropods (identity and abundance) and
 leaf damage (e.g., Karp et al. 2016). Visual arthropod surveys were
 conducted between 09:00 and 13:00 hours, always in absence of heavy rain
 (Blubaugh et al. 2018). Arthropods were counted on each plant and
 identified to the most detailed taxonomic resolution possible (e.g., lady
 beetle, aphid, parasitzed aphid [hereafter, “aphid mummy”]). We did not
 collect arthropods to avoid altering community dynamics.
 To measure leaf damage, we photographed one randomly selected
 leaf per plant per survey period. To randomly select leaves, we used the
 RAND function in Microsoft Excel to randomly assign a number between 0 and
 1 for leaf 1 (closest to the ground) through the terminal leaf. We then
 sorted the random numbers and photographed the leaf number in the row
 containing the largest random number. We repeated the selection process
 for each of the 24 plants per farm. We photographed leaves attached to the
 plant against 1-cm graph paper in the mid-growth survey to minimize plant
 damage but removed and pressed leaves overnight prior to photographing in
 the late-growth survey when the plants were generally large and well
 established. Using the images, we calculated the percent of leaf damaged
 using the software ImageJ (Kalka et al. 2008, Rueden et al. 2016). We
 estimated crop yields by cutting all mature heads 14 cm below the crown
 height following the standard of farms in the study and immediately
 measured mass in grams. Heads were harvested from 26 Aug to 21 Oct 2018,
 depending on the variety planted, grower management, and microclimate of
 experimental plots. We examined the impacts of bird
 exclusion (i.e. treatment), survey period (mid-growth vs. late-growth),
 and percent natural cover in the landscape (1000 m) on 1) crop pest
 abundance (aphids, caterpillars, and flea beetles), 2) arthropod natural
 enemy abundance (aphid mummies, all predators, spiders, and syrphid
 larvae), and 3) leaf damage using generalized linear mixed effects models
 via the <em>glmmTMB</em> package in R (Magnusson et al. 2016).
 We summed arthropods across the four plants per plot and used plot as the
 unit of replication in our arthropod analyses. We averaged leaf damage
 across the four plants per plot for leaf damage analyses. We included plot
 nested within farm as a random effect in all models. Models examining
 abundance used a negative binomial distribution, and models examining leaf
 damage used a betabinomial distribution, both to account for
 overdispersion in the data. To examine the importance of bird exclusion,
 survey period, and percent natural cover in the landscape, we compared all
 nested models from our global model that included their three-way
 interaction. We constructed a complete model set consisting of an
 intercept-only null model as well as all single variable, additive, and
 interactive models with the ‘treatment’, ‘% natural’, and ‘survey period’
 variables as fixed effects. All 15 models in our set included plot nested
 within farm as a random effect. We compared models
 based on AICc in the <em>bbmle </em>package in R (Bolker 2020)
 and identified the top-competing models as those within ∆AICc &lt; 2.0
 of the most well-supported model (Anderson et al. 2001, Burnham and
 Anderson 2002). We then estimated covariate effects by model averaging
 among the best-supported models (within 2 ∆AICc of the best-supported
 model) using the model.avg function in the <em>MuMIn
 </em>package in R (Burnham & Anderson 2002, Barton 2020). We
 considered covariates as strong predictors of the response variables if
 they appeared in the top models (∆AICc &lt; 2.0) and their
 model-averaged 95% confidence intervals did not overlap zero.
 To compare broccoli head weight (“yield”) between treatments, we
 used the log-response ratio as an effect size metric (Hedges et al. 1999).
 We averaged the head weight across plants within treatments within farms
 and calculated the effect size by taking the natural log of the average
 head weight of controls, divided by the average exclosure plot head
 weight. We used the log-response ratio rather than absolute head weight
 per plot to account for variation in crown weight between varieties grown
 by different farms. This is because varieties can differ greatly in weight
 (min = 53 g [“broccolini” variety] to max = 356 g [“gypsy” variety]). We
 installed cages within commercial fields planted by cooperating farms, so
 we were unable to mandate that all growers plant the same variety. We used
 a two-tailed one-sample t-test to determine whether the mean head weight
 differed between bird exclosure and bird access plots. We then used a
 linear regression model to assess whether difference in head weight by
 treatment varied across the landscape gradient.
 <strong>Molecular Diet Analysis</strong>
 We collected fresh fecal samples from mist-netted birds for
 molecular diet analysis. We visited farms in a south–north transect and
 generally visited each farm twice per year, but variation occurred due to
 weather and grower schedules. In total, we collected 2042 feces from 2024
 birds from 76 species captured on 47 farms between 27 Apr – 18 Sep 2016
 and 03 May – 29 Aug 2017 (Figure 1A; Figure 2D). We note that the
 exclosure experiment arthropod survey dates coincide with the “late”
 mist-net season and occurred the following year (2018). This was due to
 logistical constraints in our ability to conduct the experiment and
 netting concurrently and to maximize the number of farms able to
 participate in the exclusion experiment in the study region where the peak
 transplant period is late-June/early-July. From the
 2024 birds captured through mist-netting, we selected five species for
 molecular diet analyses (American Goldfinch, <em>Spinus
 tristis</em>; American Robin, <em>Turdus
 migratorius</em>; Savannah Sparrow, <em>Passerculus
 sandwichensis</em>; Song Sparrow, <em>Melospiza
 melodia</em>; and White-crowned Sparrow, <em>Zonotrichia
 leucophrys</em>). We selected these species because they were
 commonly observed foraging in broccoli plantings (Smith et al. 2020a;
 Smith et al. <em>in press</em>) and had sufficient numbers (≥
 50) of fecal samples to model variation across farms (Supplementary
 Material Dataset S1). Between the five focal species, we had 357 fecal
 samples from birds captured on 34 of the farms. The number of samples per
 species ranged from 53 (American Robin) to 103 (White-crowned Sparrow).
 Our choice of focal species was constrained by relatively low captures per
 bird species (see Supplementary Materials Dataset S1 for full details on
 number of samples per bird species). Thus, we included all species with
 &gt; 50 fecal samples in our diet analyses, provided that they were
 also observed foraging in brassica crops. To determine if species foraged
 in brassica crops, we conducted formal point count surveys from 2016–2019
 described in Smith et al. (2020a, <em>in press</em>). In
 addition, individual birds observed using broccoli fields during the
 experiment setup, arthropod surveys, and yield surveys were identified and
 recorded (Supplementary Material Dataset S1). Our
 mist-net protocol was previously described in detail in Smith et al.
 (2020a). Briefly, we placed 4–8 mist-nets around farms in locations
 selected to maximize capture rates and moved nets if capture rates were
 low and high activity was noted in another location. We occasionally
 placed nets immediately adjacent to each other to create longer nets along
 habitats with high bird abundances. Captured birds were placed in clean
 cloth bags, given unique leg bands, weighed, and measured. Feces were
 placed in 200 proof ethanol in cryotubes and immediately stored in a
 liquid nitrogen shipment tank until shipment to Washington State
 University on dry ice. Samples were stored at -80ºC until DNA
 extraction. To extract prey DNA from avian feces, we
 used QIAamp® DNA stool mini kits (Qiagen, Hilden, Germany) following the
 manufacturer’s protocol with modifications to increase DNA yields
 described in Supplemental Methods of Smith et al. (2020a). Each bird
 sample was then PCR amplified for arthropod specific CO1 regions using the
 ZBJ-ArtF1c (F) (AGATATTGGAAC*TTATATTTTATTTTTGG) and ZBJ-ArtR2c
 ®(*ACTAATCAATT*CCAAATCCTCC) primer set (Zeale et al. 2011). Arthropod
 primers were used with the following protocol: 95°C for 5 min, followed by
 a touchdown protocol of 16 cycles of 94°C for 30 sec, 61°C (decreasing
 0.5°C/cycle) for 30 sec and 72°C for 30 sec, then 17 cycles of 94°C for 30
 sec, 53°C for 30 sec and 72°C for 30 sec followed by a final extension at
 72°C for 10 min. PCR products were cleaned using AMPure beads following a
 modification of Rohland and Reich (2012), then pooled for each sample.
 Pooled samples were then prepared for Illumina sequencing following the
 protocol for NEBNextUltra II DNA Library Prep Kit for Illumina. After end
 repair and adaptor ligation, we again cleaned amplicons using AMPure
 beads. We then used NEBNext Multiplex Oligos for Illumina Index Primer
 Sets 1 and 2 for primer enrichment. All samples were cleaned a final time
 using AMPure beads, quantified using a Qubit Fluorometer (ThermoFisher
 Scientiﬁc) to quantify the DNA concentration of each library, then pooled
 and submitted for sequencing on the Illumina MiSeq platform (2x300bp) at
 the University of California Riverside Institute for Integrated Genome
 Biology. We included negative controls during the sequencing process to
 detect contamination or sequencing errors but did not detect any
 issues. Demultiplexed sequences were processed using
 the Anacapa Toolkit (Curd et al. 2018, Curd & Ogden 2018). First,
 sequences were trimmed using cutadapt (Martin 2011) and quality sorted
 using the fastx-toolkit (Gordon & Hannon 2010). dada2 (Callahan et
 al. 2016) was used for denoising, dereplicating, merging and removing any
 chimeric sequences. Taxonomy was assigned using Bowtie 2 (Langmead
 & Salzberg 2012) and a Bayesian Least Common Ancestor algorithm.
 For each bird, we then assessed the total number of reads for each
 Amplicon Sequence Variants (ASV) summarized at a bootstrap confidence of
 95. Because the number of reads does not always
 correlate to the amount of tissue of a specific diet item in a sample
 (Bell et al. 2019, Guenuning et al. 2019), we analyzed the percentage of
 fecal samples that contained DNA from various prey items. We restricted
 analyses to samples with at least 100 identified reads and normalized
 sequenced read counts by dividing the number of reads per amplicon
 sequence variant (ASV) by the total number of species-level identified
 ASVs assigned to phyla Arthropoda and Mollusca. We considered species to
 be present if they represented at least 1% of the reads per sample and at
 least 5 reads per ASV (Garfinkel et al. 2020). Prey items were identified
 as either brassica pests (<em>Plutella xyllostella</em>,
 <em>Trichoplusia ni</em>, <em>Pieris rapae</em>,
 <em>Phyllotrea striolata</em>, <em>Brevicoryne
 brassicae</em>, and <em>Myzus persicae</em>) or natural
 enemies observed to be present on farms in our study and/or on farms in
 other studies done in similar systems (Blubaugh et al. 2018, 2021). In
 addition, while pest sequences were identified to the species level, we
 conducted our analyses at the family level (i.e. Plutellidae, Noctuidae,
 Pieridae, Chrysomeldae, and Aphididae) to align with natural enemies. We
 had to aggregate natural enemies at the family level due to our inability
 to accurately verify species-specific presence of these taxa in our
 system. Other taxa detected (e.g., detritivores or other species not known
 to be brassica pests or arthropod natural enemies of brassica pests) were
 excluded from our analysis. Predation pressure on
 trophic groups was calculated by assigning detected prey species as either
 in the family of known natural enemies of brassica pests (predators and
 parasitoids) or in the family of known brassica pest species. We examined
 the impacts of bird species, season (late spring/early summer survey or
 late summer/early fall survey), and percent seminatural cover in the
 landscape on the likelihood of individual birds having 1) arthropod
 natural enemy taxa and 2) brassica pest taxa detected in their fecal
 samples. To do so, we used generalized linear mixed effects models with a
 binomial distribution via the <em>glmmTMB</em> package in R
 (Magnusson et al. 2016). We included farm as a random effect in all
 models. Models were run using the Bobyqa optimizer. We repeated the model
 selection approach described above for our experiment. We constructed a
 complete model set consisting of an intercept-only null model as well as
 all single variable, additive, and interactive models with the ‘species’,
 ‘% natural’, and ‘season’ variables as fixed effects. All 14 models in our
 set included farm as a random effect. For these models,
 we assessed if variables improved model fit using likelihood ratio tests
 (due to species being a categorical variable) for all models with weights
 &gt; 0.05. We used generalized Tukey HSD tests in the
 <em>multcomp </em>package in R (Hothorn et al. 2008) to
 examine differences in categorical predictor variables that had high
 support (were included in models with &lt; 2 ∆AICc) and improved model
 fit (likelihood ratio tests).