Skip to main content
Dryad logo

Natural enemy-herbivore networks along local management and landscape gradients in urban agroecosystems


Philpott, Stacy et al. (2020), Natural enemy-herbivore networks along local management and landscape gradients in urban agroecosystems, Dryad, Dataset,


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.


Study system

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 m2 to 15,400 m2. 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.

Natural enemy and herbivore surveys and identification

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 B. oleraceae plants within a 20 x 20 m plot located in the center of each garden. If there was more than one variety of B. oleraceae 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 (Brevicoryne brassicae), 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.

Vegetation, ground cover, and landscape characteristics

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 B. oleraceae plants within the 20 x 20 m plot, and the number of B. oleraceae 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 vegan 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 <5% of the surrounding landscape and were not included. We used the vegan package in R to calculate landscape diversity (e.g. modified Shannon-Wiener diversity index (H’)) for each garden at the 2 km scale. 

Natural Enemy-Herbivore Network Construction

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 B. oleraceae 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.

We examined three network metrics: estimated interaction richness, vulnerability, and trophic complementarity. We chose these three metrics because there is empirical or theoretical evidence that pest control functions should change as these metrics change. Specifically, interaction richness measures the ecological redundancy and robustness of a community. Vulnerability is a measure of the average number of natural enemy species per prey species. Trophic complementarity 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 vegan 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 bipartite 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.


Data analysis

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 B. oleracaea 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 B. oleracaea 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 car 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 MuMIn 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 visreg 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.


National Institute of Food and Agriculture, Award: 2016-67019-25185