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.

Study Sites

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.

Landscape context

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

 Bird Exclusion Experiments

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 glmmTMB 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 bbmle package in R (Bolker 2020) and identified the top-competing models as those within ∆AICc < 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 MuMIn 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 < 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.

Molecular Diet Analysis

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, Spinus tristis; American Robin, Turdus migratorius; Savannah Sparrow, Passerculus sandwichensis; Song Sparrow, Melospiza melodia; and White-crowned Sparrow, Zonotrichia leucophrys). We selected these species because they were commonly observed foraging in broccoli plantings (Smith et al. 2020a; Smith et al. in press) 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 > 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, in press). 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 Scientific) 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 (Plutella xyllostella, Trichoplusia ni, Pieris rapae, Phyllotrea striolata, Brevicoryne brassicae, and Myzus persicae) 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 glmmTMB 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 > 0.05. We used generalized Tukey HSD tests in the multcomp package in R (Hothorn et al. 2008) to examine differences in categorical predictor variables that had high support (were included in models with < 2 ∆AICc) and improved model fit (likelihood ratio tests).