Study region

The study was carried out from June to October 2015 in the Finger Lakes Region (42°26′N, 76°30′W) of New York State, USA. The landscape in this region is characterized by a mosaic of cropland and semi-natural habitats. Cropland in these landscapes mainly consisted of corn, soybean, winter wheat and crucifers, while semi-natural areas are composed of shrublands, woody wetlands, and mixed forest. We selected 11 farms across the study area to encompass a gradient of landscape complexity from landscapes with large amounts of semi-natural habitats (2% cropland) to simple landscapes dominated by crops (50% cropland) within a 1000 m radius around each farm. All farms selected for the study were either organic or used minimal inputs for pest management.

To quantify the landscape composition surrounding each farm, the proportion of semi-natural areas and cropland were calculated at three scales: 500 m, 1000 m, and 2000 m. These spatial scales are suitable for analyzing the effects of landscape context on pest control and natural enemies21. The landscape was characterized using the 2015 National Agricultural Statistics Service Cropland Data Layer for New York78 in ArcGIS 10.1.

Experimental plots

Seeds of fresh-market cabbage (B. oleracea var. capitata cv. Capture) were grown in an organic potting mix (sunshine®, Sun Gro Horticulture Inc., Bellevue, WA, USA), and fertilized with organic fish fertilizer 2-4-1(N-P-K) (Neptune’s Harvest®, Gloucester, MA, USA) for seven weeks under greenhouse conditions. Plants were eight weeks old when they were transplanted to the field.

On each of the 11 farms, we established two experimental plots. One plot was randomly chosen for the augmentative predator release treatment while the other served as a non-release control. Plots within the same farm were separated by 334 ± 41 m (mean ± 1 SE), and the mean distance between farms was 7.2 ± 2.3 km. Care was taken to minimize fine-scale landscape heterogeneity between experimental plots within the same farm. Plots within the same farm primarily differed in the predator release treatment, while landscape context, plot size and shape, and abiotic conditions were similar for each pair.

Each experimental plot consisted of ten 7.2-m rows, with 15 cabbage plants per row. Row and plant spacing were 0.9 m and 0.45 m, respectively. Plants were transplanted across farms over two consecutive weeks in mid-June 2015. Experimental fields within the same farm were planted on the same day. Plants were fertilized during transplanting and again one month later using 8-3-3 (N-P-K) granular compost at a rate of 5 kg/100 m2 (Kreher’s® composted poultry manure, Clarence, NY, USA). All experimental plots were managed without fungicides or insecticides, and weeds were removed at two-week intervals.

Augmentative releases of predators

The predator release treatment included both Podisus maculiventris nymphs and Hippodamia. convergens adults. Both the nymphal and adult stinkbugs display high predation rates on lepidopteran larvae, so we released fourth and fifth instars in our experiments to minimize dispersal after release and increase the potential for season-long pest control. Ladybird beetle larvae were not available commercially, which precluded us from using less-mobile stages. Predators were released three times throughout the season at the seedling, pre-cupping, and early head formation growth stages79. Releases were conducted early in the season, as previous studies have shown that early control is key to the success of biocontrol strategies in field settings80,81. Approximately 200 stinkbugs and 600 ladybeetles were released per plot each time by carefully deploying them on the leaves. These release rates equaled 1.3 sting bug nymph/plant and 4 ladybird adults/plant. These are commonly recommended release rates by commercial vendors82,83,84,85. No predators were released in the control plots.

P. maculiventris were obtained from eggs purchased from a commercial supplier (Beneficial Insectary Inc., Redding, CA, USA) and reared on a diet of mealworms, Tenebrio molitor (L.), and cabbage seedlings. The stinkbug colony was kept at 25.5 ± 2.0 °C, 60% RH, and a photoperiod of 16:8 (L:D) following the methods of De Clercq et al.86. Adults of H. convergens were obtained from a commercial supplier (Arbico Organics, Oro Valley, AZ, USA). Ladybird beetles were stored at 7 °C until we released them in the field.

Measuring predation rates in the field

Predation rates provided by resident and augmented predators were quantified using Trichoplusia ni larvae and eggs as sentinel prey. Trichopluisia ni were commercially available and easier to manipulate in field studies than P. rapae. Trichopluia ni larvae and eggs were obtained from a commercial insectary (Benzon Research Inc., Carlisle, PA, USA). For estimating larval predation, 5 third-instars (13.6 ± 0.23 mm long) were placed on the upper leaves of four randomly selected plants per plot (i.e., 20 larvae per plot). After 24 h of exposure in the field, the remaining larvae were counted to determine the number of larvae consumed by predators. Larvae were considered predated if they were completely missing, or showed evidence of predation such as necrotic tissue around an open wound.

To estimate egg predation per plot, paper discs containing approximately 30 T.ni eggs (range: 19–76) were fixed to the underside of 10 × 10 cm pieces of corrugated plastic board (Coroplast®, Vanceburg, KY, USA) that provided a standardized foraging platform for predators. Five egg platforms were positioned at crop height and placed between the leaves of the plants where sentinel larvae were deployed. All egg masses were inspected after 24 h, and the number of eggs remaining were counted to determine predation rates. Eggs were considered predated if they were missing, the chorion presented clear evidence of attack (i.e., chewing predator), or the contents of the egg had vanished (i.e., attacks by a piercing-sucking predator). To distinguish sentinel prey predation from unknown losses due to handling and rainfall, we enclosed one plant per site in a cage that excluded natural enemies. Cages consisted of a 0.2 × 0.2 × 0.2 m3 mesh plastic screen (BioQuip, Rancho Dominguez, CA, USA) with openings of 1.1 × 0.7 mm2, and whose bottom edges were buried 5 cm into the ground. Plants in these cages were infested with sentinel prey in the same fashion as the uncaged plants. Net mortality due to predation was determined by assessing mortality from uncaged plants and subtracting it from mortality from caged plants.

We repeated the sentinel prey experiment three times per plot at the seedling, pre-cupping, and early head formation growth stages. Thus, we had three temporally separated dates that allowed us to account for the temporal differences in predation rates throughout the season.

Sampling of lepidopteran pests and their natural enemies

To assess lepidopteran abundance, plants were visually inspected for larvae during the seedling, pre-cupping, early head formation, and maturation growth stages79. In each plot, ten randomly selected plants were destructively sampled and the number of larvae were recorded on each plant. To avoid possible edge effects, plants within 1 m of the edge of the plot were not sampled. A total of 294 caterpillars were collected in the experimental plots, with P. rapae as the dominant species (94% of the total caterpillars collected) followed by P. xyllostela (5%) and T. ni (0.4%).

Naturally occurring predators and parasitoids were sampled using yellow sticky cards, pitfall traps, and visually inspected plants. Natural enemies from these samples were categorized into functional groups as foliar-foraging predators, ground-dwelling predators, and parasitoids. Our analysis was restricted to species known to attack either lepidopteran eggs or larvae based on previous observational and experimental studies (e.g.43,87). For foliar-foraging predators, we focused on the three dominant species of coccinellids in our system: the native Coleomegilla maculata and the two-exotic species, Harmonia axyridis and Propylea quatuordecimpunctata. Abundance of all three coccinellid species were pooled to obtain the overall abundance of relevant foliar-foraging predators for each plot. For ground-dwelling predators, carabid beetles were collected and identified to species. Following identification, we gathered information from the literature to further classify carabids into three diet categories: carnivorous, omnivorous or phytophagous88,89. Only carnivorous species were kept in further analyses. Altogether 25 predatory carabid species were collected, of which three species (Bembidion quadrimaculatum, Poecilus chalcites and Poecilus lucublandus) made up 66% of the total capture. As with coccinellids, the abundance of predatory carabids for each plot was pooled in subsequent analysis. Lastly, we measured parasitoid abundance by focusing our sample efforts on Cotesia rubecula (Hymenoptera: Braconidae), the most important specialist parasitoid of P. rapae larvae in the study region90. The parasitoids of the T. ni and P. xylostella were not investigated because both pests occurred in small numbers in our system (i.e. <6% of the total caterpillars collected).

Sampling for all natural enemies was conducted three times during the season at the seedling, pre-cupping, and early head formation stages. Foliar-foraging predators, ground-dwelling predators and parasitoids were sampled using sticky cards, pitfall traps and visual inspection of plants, respectively. On each sampling time, one sticky card (15 × 30 cm, BioQuip, Rancho Dominguez, CA, USA) was positioned at crop height in the center of each plot. The sticky cards were retrieved after 15 days and the number of foliar-foraging predators were recorded. For the pitfall traps, a 540 mL clear plastic cup (9 cm diameter openings, Fabri-kal corp., Kalamazoo, MI, USA), was filled with a mixture of water and a few drops of organic, odorless detergent (Dr. Bronner’s Unscented Pure Castile Soap, Vista, CA, USA). A total of five traps were placed within the rows between cabbage plants; four traps were located near the corners and one in the center row of the plot. Each trap was protected from rain and direct sunlight by a plastic plate (15 cm in diameter) held approximately 10 cm above the trap. Pitfall traps were collected after 24 h and the number of ground-dwelling predators was recorded. Finally, on each sampling date, parasitoid abundance was estimated by counting the total number of parasitoid cocoons (i.e, pupa) on ten randomly selected plants per experimental plot. Parasitoids were identified using diagnostic morphological characters described by Van Driesche (2008)91.

Plant damage and crop biomass

Insect damage and crop biomass were assessed from the same ten plants used for lepidopteran censuses at four sampling times during the season. Damage was quantified using a modified version of the method of Lim et al.92, where a plant is classified into one of the following eight categories based on the percentage of leaf area removed: 0, 5–10,10–20, 20–40, 40–60, 60–80, 80–100, or 100%). Visual estimates of damage provide the fastest and most cost-effective method for quantifying herbivory93, and previous studies have shown they can be precise and accurate to estimate economic thresholds for lepidopteran defoliation in cabbage94. For analysis, we assumed the estimated proportion of damage on each plant to be the median of each category (0, 7, 15, 30, 50, 70, 90, 100, respectively). Crop biomass was determined by weighing the plants after they had been oven-dried at 60 °C for 7 days. Although the crop biomass at the end of the season (i.e. maturation growth stage) is a measure of crop yield (i.e. marketable cabbage head weight), we used the crop biomass throughout the season in our analyses rather than only final biomass, as the former allowed us to account for the temporal differences in the effectiveness of augmentative biocontrol. Analysis using only final crop yield produced qualitative similar results (Supplementary Fig. S1).

Laboratory experiment

Controlled lab experiments were conducted to quantify the individual and combined effect of stinkbugs and ladybird beetles on pest predation, independently from the effects of landscape context. Experimental units were 28 × 28 × 28 cm cages, covered on all sides, expect the bottom, with a mesh screen opening of 1.1 × 0.7 mm (BioQuip, Rancho Dominguez, CA, USA). Each experimental unit consisted of a single potted cabbage plant (B. oleracea var. capitata cv. Capture) with six fully-expanded true leaves. To begin the experiment, all cages received 5 third-instar larvae and one egg mass (approximately 30 eggs) in the same fashion as the sentinel field experiment. T. ni larvae were allowed to establish for 1 h before predator introduction.

Predators were released into individual cages according to four treatments: control (no predators added), stinkbug treatment (2 fifth-instars added), ladybird beetle treatment (5 adults added), and the interaction treatment (2 fifth-instar stinkbugs and 5 ladybird beetles adults added). These densities were chosen because they mimicked those used in the sentinel field experiment. Each treatment was replicated eight times. The experiment had an additive design (i.e., overall predator density is higher in the multi-species treatment compared to the single-predator treatment), because this approach better reflects the effect of predator augmentation. Predators were starved for 24 h before being introduced to standardize hunger levels across treatments and then allowed to feed for an additional 24 h, after which the number of larvae and eggs remaining in the cages were recorded. The experiments were conducted at 25 ± 2 ◦C, 60 ± 5% RH and a 14:10 (L:D) h photoperiod.

Statistical analysis

To analyze the direct effect of the abundance of naturally occurring enemies on biocontrol of sentinel prey, pest incidence, and plant damage, we used linear mixed-effect models in R with the nlme package95. Abundances were averaged for each functional group separately (i.e. foliar-foraging predators, ground-dwelling predators, and parasitoids) and for each sampling period. Response variables were square-root-transformed to meet assumptions of normality and homoscedasticity96. For all models, we also included farm as random effect to account for other potential sources of variability associated with each geographic location (e.g., environmental or management intensity differences). Statistical significance of the abundance of each functional group was assessed by conditional F-tests97.

The effects of landscape complexity and potential interactions with predator releases, on lepidopteran pest abundance, natural enemy abundance, predation rates, plant damage, and crop biomass, were examined using linear (lmer) and generalized linear mixed-effect models (glmer)98. Fixed factors in the models included treatment (with or without predator releases), landscape complexity, and the treatment by landscape interaction. Landscape complexity was defined as either the proportion of cropland or the proportion of semi-natural areas as both variables were highly correlated at all scales (Spearman’s r s < −0.45, P < 0.001 at all scales). Random effects in all models included farm and sampling time to account for the crossed experimental design (i.e., each plot was measured on multiple dates and multiple plots were measured on each date). Response variables were square-root transformed prior to analysis to meet normality assumptions and avoid heteroscedasticity. Assumptions were checked according to the graphical validation procedures recommended by Zuur et al.96. Models for foliar-foraging predators did not meet distributional assumptions, and therefore were analyzed using a generalized linear mixed model with a Poisson error distribution. Model simplification was done using a backwards-stepwise selection (lmer) or likelihood ratio tests (glmer) based on Akaike’s Information Criterion (AIC) where non-significant predictors were removed (P > 0.05). We assessed the statistical significance of fixed effects and interaction terms by F-tests based on the Satterthwaite approximation (lmer) or Wald Z-test (glmer)99,100. Separate models were fitted for each landscape scale (i.e., 500, 1000, and 2000 m), and the scale with the highest explanatory power for each response variable was determined by comparing the AIC values of the minimum adequate models101. The most predictive scale for each response variable was then used in further analyses. Subsets of best models for each response variable are provided in Supplementary Table S1. Mantel test102 indicated no spatial autocorrelation in the residuals of the final models (Supplementary Table S2).

To better understand potential differences between the predator release treatment and the control, we used the final models to estimate the marginal means and 95% confidence intervals for each response variable with the “emmeans” package in R103. For all response variables (i.e., lepidopteran abundance, predation rates, plant damage, crop biomass, and natural enemy abundance), we used preplanned contrast to determine whether mean differences between plots with and without predator releases (i.e, effect size) were significant. We first estimated the mean effect size across the entire landscape complexity gradient to get an overall quantitative assessment of the consequences of predator augmentation for each response variable. In a second group of comparisons, we estimated the effect size of each response at even intervals over the landscape complexity gradient (range: 0–0.6) to test the hypothesis that the effects of augmented predators were contingent on the characteristics of the surrounding landscape. Pairwise multiple comparisons were calculated using the Bonferroni correction for an overall error rate of 0.05. Comparisons were conducted using the emmeans package.

For the laboratory experiment, we examined predation rates on lepidopteran larvae and eggs using a two-way ANOVA followed by a Tukey HSD test at P < 0.05 including the factors: stinkbugs (with or without), ladybird beetles (with or without), and their interaction. Predation rates were log-transformed to meet the assumptions of the analysis. To further examine these data, we used a multiplicative risk model104 followed by ANOVA comparing the expected and actual predation rate values to assess whether combined predators act independently (i.e., observed and predicted values do not differ, so that its combined effect is additive), antagonistically (i.e., observed values are less than the predicted values), or synergistically (i.e., observed values exceed the predicted values) on prey populations60. All statistical analyses were done using R v. 3.2.3105.