DOI: 10.7291/d1dh3n
title: Natural enemy-herbivore networks along local management and landscape gradients in urban agroecosystems

description: Ecological networks can provide insight into how biodiversity loss and
 changes in species interactions impact the delivery of ecosystem services.
 In agroecosystems that vary in management practices, quantifying changes
 in ecological network structure across gradients of local and landscape
 composition can inform both the ecology and function of productive
 agroecosystems. In this study, we examined natural enemy-herbivore
 co-occurrence networks associated with Brassica oleracea (cole crops), a
 common crop in urban agricultural systems. Specifically, we investigated
 how local management characteristics of urban community gardens and the
 landscape composition around them affect a) the abundance of B. oleracea
 herbivores and their natural enemies, b) the natural enemy: herbivore
 ratio, and c) natural enemy-herbivore co-occurrence network metrics. We
 sampled herbivores and natural enemies in B. oleracea plants in 24
 vegetable gardens in the California central coast region. We also
 collected information on garden characteristics and land-use cover of the
 surrounding landscape (2 km radius). We found that increased floral
 richness and B. oleracea abundance were associated with increased
 parasitoid abundance, non-aphid herbivore abundance, and increased network
 vulnerability; increased vegetation complexity suppressed parasitoid
 abundance, but still boosted network vulnerability. High agricultural
 land-use cover in the landscape surrounding urban gardens was associated
 with lower predator, parasitoid, and non-aphid herbivore abundance, lower
 natural enemy: herbivore ratios, lower interaction richness, and higher
 trophic complementarity. While we did not directly measure pest control,
 higher interaction richness, higher vulnerability, and lower trophic
 complementarity are associated with higher pest control services in other
 agroecosystems. Thus, if gardens function similarly to other
 agroecosystems, our results indicate that increasing vegetation
 complexity, including trees, shrubs, and plant richness, especially within
 gardens located in intensively farmed landscapes, could potentially
 enhance the biodiversity and abundance of natural enemies, supporting
 ecological networks associated with higher pest control services.
description: <b><i>Study system </i></b>
 We studied natural enemy (predator and parasitoid) and herbivore
 abundance and co-occurrence networks in urban gardens in the California
 central coast region between May and August 2017. We selected 24 urban
 community gardens in Monterey (7 gardens), Santa Cruz (9 gardens), and
 Santa Clara (8 gardens) counties for field research (Fig. 1). The sites
 within these three counties are distributed in two different California
 ecoregions - the Monterey Bay Plains and Terraces (Santa Cruz and Monterey
 County sites) and the Bay Terraces/Lower Santa Clara Valley (Santa Clara
 County sites), which are broadly different in terms of biotic and abiotic
 phenomena that may influence the ecology of those systems. All gardens are
 community gardens (with vegetables, fruit trees, and ornamental plants)
 managed collectively or in individual allotments (plots), and the gardens
 range in size from 444 m<sup>2</sup> to 15,400
 m<sup>2</sup>. Each garden had been in production for between
 2-50 years during the time of the study. All gardens were separated from
 each other by a minimum of 2 km.
 <b><i>Natural enemy and herbivore surveys and
 identification</i></b> We surveyed gardens
 four times over the summer growing season for natural enemies and
 herbivores (between 15-19 May, 19-23 June, 17-20 July, and 14-17 August
 2017). We haphazardly selected up to 20 <i>B. oleraceae</i>
 plants within a 20 x 20 m plot located in the center of each garden. If
 there was more than one variety of <i>B. oleraceae</i> in the
 plot (e.g., broccoli, curly kale, cabbage), or more than 20 individual
 plants, we selected plants to represent the relative abundance of
 different varieties and spatial distribution across the plot. We measured
 the height and width, the flowering status , and the variety (e.g. curly
 kale, broccoli) of all surveyed plants. We visually examined all above
 ground plant parts (e.g. leaves, stems, fruits, flowers) of each plant for
 arthropods. We then manually collected all arthropods encountered either
 directly on the plant or hovering above the plant using forceps and vials.
 The exception was cabbage aphids (<i>Brevicoryne
 brassicae</i>), which we counted but did not collect due to their
 high densities and easy identification. All collected arthropods were
 placed into a vial with 70% ethanol and individually marked.
 We identified all arthropods to order, family, genus, or species
 as necessary to determine the trophic guild of the arthropod (e.g.,
 herbivore, parasitoid, predator). We identified the trophic guilds. In
 cases where we failed to collect an arthropod (e.g. it flew away before
 collection), we described the arthropod and noted the order or family as
 possible. Only those arthropods that were collected or described at least
 to order were included in the analyses.
 <b><i>Vegetation, ground cover, and landscape
 characteristics</i></b> Within the 20 x 20
 m plot at the center of each garden, we sampled vegetation and ground
 cover characteristics. We sampled canopy cover with a concave spherical
 densiometer at the center of each plot, and 10 m to the N, S, E, and W of
 the center. We counted the number of trees and shrubs, the number of tree
 and shrub species, the number of trees and shrubs in flower, and the
 number of tree and shrub species in flower. We also counted the total
 number of <i>B. oleraceae</i> plants within the 20 x 20 m
 plot, and the number of <i>B. oleraceae</i> plants sampled
 during each visit. Within the 20 x 20 m plots, we randomly selected eight
 1 x 1 m plots within which we identified all herbaceous plants (except
 grass) to morphospecies; counted the number of flowers and number of
 species in flower; measured the height of the tallest herbaceous
 vegetation; and visually estimated the percent of the plot covered by a)
 bare ground, b) grass, c) herbaceous plants, d) rocks, e) leaf litter, f)
 straw, and g) mulch or wood chips. We took vegetation and ground cover
 data on the same days that we sampled arthropods, and all values for each
 site (except for herbaceous species richness) were averaged across the
 four sample dates. For herbaceous species richness, we estimated the total
 herbaceous plant species richness (Chao1) in each site across all sample
 dates with the estimateR function in the <i>vegan</i> package
 for R. We calculated a vegetation complexity index (VCI) for each site. To
 calculate the VCI, we included canopy cover, number of trees and shrubs,
 number of tree and shrub species, the number of trees and shrubs in
 flower, the number of trees and shrubs and tree and shrub species in
 flower, estimated species richness of herbaceous plant species, herbaceous
 plant cover, and height of the tallest vegetation. We scaled values for
 each variable from 0 to 1 by dividing by the highest value measured across
 all sites. We then averaged values for the eight variables to yield an
 overall VCI between 0 (low vegetation complexity) and 1 (high vegetation
 complexity). Thus, overall, we collected or calculated a total of 18
 vegetation and ground cover variables in each garden.
 We used land-cover data from the 2011 National Land Cover
 Database (NLCD, 30 m resolution) and calculated the percentage of
 land-cover types in 2 km buffers from the center of each garden. We
 created four land-cover categories: 1) natural (including the NLCD
 categories of deciduous, evergreen, and mixed forests, dwarf scrub,
 shrub/scrub, and grassland/herbaceous), 2) open (including lawn grass,
 park, and golf courses), 3) urban (including low, medium, and high
 intensity developed land), and 4) agriculture (including pasture/hay and
 cultivated crops). Other land-cover types covered &lt;5% of the
 surrounding landscape and were not included. We used the
 <i>vegan</i> package in R to calculate landscape diversity
 (e.g. modified Shannon-Wiener diversity index (H’)) for each garden at the
 2 km scale.  <b><i>Natural Enemy-Herbivore
 Network Construction</i></b> We constructed
 natural enemy networks based on positive co-occurrence of families of
 natural enemies (predators or parasitoids) and herbivores on the same
 plant (as in Bell et al. 2010). For each site, we calculated the mean
 number of aphids, non-aphid herbivores, predators, parasitoids, and
 natural enemies (predators plus parasitoids) per <i>B.
 oleraceae</i> plant. In order to estimate feeding links between
 individual families of natural enemies and herbivores, we used qualitative
 diet and host information for all the observed predator and parasitoid
 taxa from the literature. If we could not find reports of interactions
 between two families, and determined it unlikely that an interaction was
 not possible due to mouth-part or physiological constraints, we marked
 this as a ‘forbidden’ link and used a different notation in our dataset to
 build the network and to differentiate from true zeroes.
 <br> We examined three network metrics: <i>estimated
 interaction richness</i>, <i>vulnerability</i>, and
 <i>trophic complementarity</i>. We chose these three metrics
 because there is empirical or theoretical evidence that pest control
 functions should change as these metrics change. Specifically,
 <i>interaction richness</i> measures the ecological redundancy
 and robustness of a community. <i>Vulnerability </i>is a
 measure of the average number of natural enemy species per prey species.
 <i>Trophic
 complementarity</i> represents the degree to which
 natural enemies share prey resources, and is often used as a predictor of
 pest control.
 To calculate estimated interaction richness, we first noted
 each positive co-occurrence of natural enemy-herbivore family pairs in
 each site, tallied abundance as the number of plants on which that family
 pair co-occurred, and then used the estimateR function with the
 <i>vegan</i> package in R to calculate Chao1 (abundance-based)
 estimated interaction richness for each site. Use of estimated interaction
 richness, rather than simple numbers of interactions, buffers against
 lower sample sizes in network studies. We used the
 <i>bipartite</i> package in R to calculate vulnerability, a
 measure of the mean number
 of natural enemy families per herbivore family and NDOF (nestedness
 metric based on overlap and decreasing fill). We then calculated
 trophic complementarity (C, a measure of the shared natural enemy families
 for each herbivore) (e.g., the inverse of NDOF where C = (100 - NODF) /
 100). For calculating network metrics, we used site as the replicate and
 the number of plants within a site on which a natural enemy and herbivore
 family co-occurred as our metric of abundance.
   <b><i>Data
 analysis</i></b> Because of the large
 number of predictor variables collected, and the potential for collinear
 variables, we first ran Pearson’s correlations to identify correlated
 variables and to select variables for subsequent analysis. We grouped
 variables into biologically relevant groups (e.g., vegetation variables,
 ground cover variables, landscape variables). We then identified variables
 that were significantly correlated with one another and selected one
 variable per group to include in further analysis. Of the correlated
 variables, we selected the one with the highest average correlation
 coefficients with other variables. We also selected variables that were
 not correlated with any others. Based on the Pearson’s correlations, we
 selected a total of nine explanatory variables: VCI, number of <i>B.
 oleracaea</i> per plot, ecoregion, garden size, percent bare ground,
 percent grass, number of flower species, agriculture within 2 km, and
 urban land cover within 2 km. We used natural log (garden size, number of
 <i>B. oleracaea</i> plants) or square root (grass cover in 1 x
 1 m plots, agriculture cover within 2 km, and urban land cover within 2
 km) transformed data for some variables to improve model fit. Because
 correlation coefficients only consider pairwise comparisons, and to assure
 that we did not have collinearity between some of the remaining
 explanatory variables, we checked the variable inflation factor (VIF) with
 the ‘vif’ function in the <i>car</i> package version 3.0-2.
 Initial tests indicated that urban cover within 2 km was collinear with
 other variables: after removing that variable from our global model, all
 VIF scores were below 2.3.  
 To examine which local and landscape factors drive
 changes in natural enemy and herbivore abundance, natural enemy: herbivore
 ratio, and network metrics, we used generalized linear models (GLMs) with
 the glm function in R. For herbivore and natural enemy variables, we
 examined aphid abundance, non-aphid herbivore abundance, predator
 abundance, parasitoid abundance, and the ratio of natural enemies to
 herbivores as dependent variables in our models. To examine changes in
 network metrics, we used estimated interaction richness, vulnerability,
 and trophic complementarity as dependent variables. For each of the
 dependent variables, we tested all combinations of the eight selected
 explanatory variables with the ‘glmulti’ function and selected the top
 model based on the AICc values. For all models where the AICc value was
 within two points of the next best model, we averaged models with the
 model.avg function in the <i>MuMIn</i> package and report
 conditional averages for significant model factors. For all dependent
 variables (except trophic complementarity), we used a negative binomial
 distribution as this provided the best fit (e.g., relatively equal
 residual deviance and df values and non-significant asymptotic Chi-square
 tests for goodness of fit). For trophic complementarity, a Gaussian
 distribution provided the best fit. We visualized all significant local
 and landscape predictors of prey removal from either top models or the
 averaged top model with the <i>visreg</i> package in R
 (Breheny and Burchett 2013). All statistical analyses were
 conducted in R version 1.1.456.
 Due to the importance of agricultural cover in the landscape in
 many models (see results), and a relatively small number of sites with
 agriculture in the landscape (n=6), we performed a subsequent analysis as
 described above, only including sites with no agriculture in the landscape
 (n=18) to determine which local and landscape factors drove variation in
 natural enemy and herbivore abundance and network metrics in those sites.
 For those models, we included VCI, number of brassicas per plot, percent
 bare ground, percent grass, number of flowering species, and urban land
 cover within 2 km to maintain VIF scores below 1.65. We did not include
 ecoregion or garden size in these models due to high
 collinearity.