Effects of Prey Distribution and Heterospecific Interactions on the Functional Response of Harmonia axyridis and Aphidius gifuensis to Myzus persicae

Natural enemy guilds normally forage for prey that is patchily distributed simultaneously. Previous studies have investigated the influence of conspecific interactions and prey distribution on the functional response of natural enemies. However, little is known about how prey distribution and heterospecific interactions between natural enemies could affect their foraging efficiency. We examined the effects of prey distribution (aggregate and uniform) and heterospecific interactions on the functional response of a predator, Harmonia axyridis (Pallas) (Coleoptera: Coccinellidae) and a parasitoid, Aphidius gifuensis Ashmead (Hymenoptera: Braconidae) to the green peach aphid, Myzus persicae (Sulzer) (Hemiptera: Aphididae). Type II functional responses were observed in all experiments. Functional response curves of single H. axyridis or A. gifuensis were higher in the aggregate treatment than in the uniform treatment when aphid densities were between 40–180 or 70–170, respectively. When comparing between aggregate and uniform treatments with the heterospecific enemy occurrence, no differences were found in the parasitism efficiency of A. gifuensis, while H. axyridis consumed more aphids in the aggregate treatment than in the uniform treatment when aphid densities were between 50–230. The functional response of individual H. axyridis was not affected by A. gifuensis under two aphid distributions. However, the functional response of a single A. gifuensis and the treatment when A. gifuensis concurrently with H. axyridis overlapped in uniform treatment of above approximately 150 aphids. Our results indicate that the predation rate of H. axyridis was affected by aphid distribution, but was not affected by heterospecific interactions. The parasitism rate of A. gifuensis was affected by aphid distribution, and by heterospecific interactions in both the aggregate and uniform treatments. Thus, to optimize the management efficiency of M. persicae, the combined use of H. axyridis and A. gifuensis should be considered when M. persicae is nearly uniformly distributed under relatively high density.


Introduction
Biological control by natural enemies is an environment-friendly and effective approach in regulating pest population, and it has received increasing research interest and has long been applied as part of integrated pest management (IPM) strategies [1][2][3]. In most ecosystems, pest species are often associated with multiple natural enemies, and more and more biological control programs have used more than one species of natural enemies [4][5][6]. However, the effect of multiple enemies in a framework to understand population dynamics of each natural enemy species and guide strategies to increase the efficacy of combining these two natural agents in biological control programs.

Plants and Insects
Chili pepper plants (Capsicum annuum L., var. 'Shulahuojian F1 , six-week-old and around 12 cm in height) were used for rearing aphids or preparing for the experiments. Myzus persicae, A. gifuensis, and H. axyridis were originally collected from chili pepper and cabbage fields at the Experimental Farm (108 • 04 18" E, 34 • 17 52" N), Northwest A&F University (Yangling, Shaanxi, China) in July 2014. Myzus persicae were maintained on chili pepper plants; A. gifuensis were originally obtained from M. persicae and cultured for at least 12 generations on M. persicae on chili pepper plants. Harmonia axyridis males and females were paired in Petri dishes (3 cm in diameter) and fed with M. persicae to allow mating and oviposition. Newly hatched H. axyridis larvae were reared individually in 3-cm Petri dishes and provided with an excess of M. persicae daily until they reached the pupal stage. Harmonia axyridis adults that emerged within 24 h were isolated and the naïve 2-3-days old unmated female adults were used in all subsequent experiments. All insect colonies and the experiments were maintained in a controlled insectary at 25 ± 1 • C, 65 ± 5% RH and a 16:8 h (L:D) photoperiod.

Functional Response of Single Predators or Parasitoids
For each predator or parasitoid species, we first measured the predation/parasitism rate of a single adult female consuming or parasitizing M. persicae on chili pepper plants. Third instar aphid nymphs were used in all experiments to avoid aphid nymph production by adults. To create third instar aphid nymphs, a pepper leaf disc (3 cm in diameter) was placed on the bottom of a small Petri dish (3 cm in diameter) with 1% agar gel. Twenty M. persicae adults were introduced in each Petri dish and removed after 24 h. The newborn aphid nymphs were maintained in the Petri dish for another 24 h. Newly molted second-instar nymphs were transferred to new leaf disks and all younger nymphs and the ecdyses were discarded. Aphid nymphs were kept and were used when they grew to the third instar.
Four chili pepper plants were arranged randomly to each ventilated cage (30 × 30 × 30 cm). Each experimental cage received 4,8,16,32,64,128, or 256 third-instar aphid nymphs. Aphids were introduced at the bottom of the stem of the chili pepper plant and allowed to acclimate for 1 h. For the aggregation treatment, all aphids were randomly allocated on one of these plants. However, for the uniform treatment, aphids were divided evenly among four plants.
Prior to the experiment, unmated adult female H. axyridis were individually transferred from the stock culture into Petri dishes (3 cm in diameter) for 24 h to standardize their hunger level. During this time, a water-saturated cotton ball was placed in each Petri dish to provide moisture. For the parasitoid, A. gifuensis mummies were collected from plants with a fine camel hair brush and placed in plastic cylindrical cages (12 cm in height by 7 cm in diameter) with 10% honey solution and inspected at regular intervals. All male and female parasitoids that emerged on the same day were placed in new plastic cylindrical cages for 24 h, and the parasitoids were left undisturbed to ensure female mating. Generally, adult males and females normally mated a few hours after emergence [37]. Mated females were used in the experiments.
Then, H. axyridis or A. gifuensis were placed individually in the center of each cage. After 24 h, predators or parasitoids were removed from the cage. For the predator treatment, the number of aphids consumed was recorded. For the parasitoid treatment, the number of parasitoid mummies was recorded after 10 days. The experiment was replicated 10 times for each treatment.

Functional Response of Paired Heterospecific Enemies
Starved H. axyridis adults and A. gifuensis mated females were prepared using the same procedures as described in the relevant sections, and the aphid density and distribution used were the same as described above. One H. axyridis adult and one A. gifuensis female were introduced to each cage. After 24 h, predators and parasitoids were removed from the experimental cages and the number of aphids preyed upon by predators was recorded. Moreover, the number of parasitoid mummies was recorded after 10 days. The experiment was replicated 10 times for each treatment.
For both the two experiments above, aphids were not replaced during the experiment. In addition, under two aphid distributions, a control treatment without predators and parasitoids was conducted with five replications for each aphid density (4, 8, 16, 32, 64, 128, or 256) to assess natural mortality rates by counting the dead aphids.

Data Analysis
All analyses were conducted using the statistical software R [47]. To evaluate the type of functional response that best fitted the data in the different experiments, a model selection and hypothesis testing was used [48]. For model selection, a logistic regression of the number of prey killed was used to identify the type of functional responses fitted with the maximum likelihood (ML) procedure. Significant negative or positive linear coefficients from the regression suggest type II or III responses, respectively [48]. When a significant negative linear coefficient from logistic regression was found, the data were then fitted to a type II functional response curve with ML estimation using the random predator Equation (1) [49], which allows for prey depletion: where N e is the number of prey eaten, N 0 is the initial prey density, a is the attack rate, T h is the handling time, and T is the total experimental duration (24 h). To compare functional response fits between natural enemies, the functional response fits were non-parametrically bootstrapped (n = 2000) to generate 95% confidence intervals (CIs) around functional response curves and the associated parameters. Equation (1) was then fitted to the bootstrapped dataset with initial parameter values that were estimated from the original ML estimates. The overlap between confidence intervals indicates that the functional responses and/or the corresponding parameters were not significantly different. Analysis of the observed functional responses modeling was carried out with the 'frair' package [50]. Data from trials of single H. axyridis, single A. gifuensis, and individual H. axyridis or A. gifuensis in heterospecific combination were analyzed with a generalized linear mixed model (GLMM) (glmer function in the lme4 package) with a binomial distribution. The dependent variables were the number of aphids killed, and the explanatory variables were aphid density and their distributions. Natural enemies tested in each replicate was treated as an observation-level random effect.

Results
In control treatments, aphid survival in both types of distribution treatment exceeded 98.5% in ventilated cages, and thus aphid's natural mortality did not attribute to background mortality.

Functional Response of Single Predators or Parasitoids
Significant negative linear terms were detected from logistic regressions for both treatments. This indicated a type II functional response for single H. axyridis or A. gifuensis ( Table 1). The attack rates and handling times of the functional response models were all significant (Table 1).
For single H. axyridis treatment, we found that aphid density, aphid distribution, and the interaction between aphid density and aphid distribution had a significant effect on aphid consumption by H. axyridis adults ( Table 2). Functional response curves overlapped at aphid densities below 40 and above 180 between aggregate and uniform treatment (Figure 1a). For single A. gifuensis treatment, the aphid density, Insects 2020, 11, 325 5 of 12 aphid distribution, and the interaction between aphid density and aphid distribution significantly affected the number of aphids parasitized by female A. gifuensis (Table 2). Functional responses of aggregate and uniform treatments overlapped at aphid densities below 70 and above 170 (Figure 1b). Bold letters indicate significant differences between treatments (p < 0.05).

Functional Response of Paired Heterospecific Enemies
Both H. axyridis and A. gifuensis exhibited type II functional responses when heterospecific enemy species were present ( Table 1). The attack rates and handling times of the functional response models were all significant ( Table 1).
The number of aphids consumed by H. axyridis was affected by aphid density, aphid distribution, and the interaction between aphid density and aphid distribution when A. gifuensis was

Functional Response of Paired Heterospecific Enemies
Both H. axyridis and A. gifuensis exhibited type II functional responses when heterospecific enemy species were present ( Table 1). The attack rates and handling times of the functional response models were all significant ( Table 1).
The number of aphids consumed by H. axyridis was affected by aphid density, aphid distribution, and the interaction between aphid density and aphid distribution when A. gifuensis was present ( Table 2). In the heterospecific enemy combination, more aphids were consumed by H. axyridis in the aggregate treatment than in the uniform treatment when aphid density was between 50 and 230 aphids ( Figure 2a). As for A. gifuensis, there was a significant effect of aphid density on the number of aphids parasitized by female A. gifuensis when H. axyridis was present (Table 2). When H. axyridis occurred, functional response of A. gifuensis overlapped across all prey densities between aggregate and uniform treatment (Figure 2b).
Insects 2020, 11, x FOR PEER REVIEW 2 of 13 present ( Table 2). In the heterospecific enemy combination, more aphids were consumed by H. axyridis in the aggregate treatment than in the uniform treatment when aphid density was between 50 and 230 aphids ( Figure 2a). As for A. gifuensis, there was a significant effect of aphid density on the number of aphids parasitized by female A. gifuensis when H. axyridis was present (Table 2). When H. axyridis occurred, functional response of A. gifuensis overlapped across all prey densities between aggregate and uniform treatment (Figure 2b). Functional response curves were overlapped between the treatment where H. axyridis was alone and the treatment where H. axyridis was sharing the experimental patch with A. gifuensis for aggregate and uniform treatments, respectively (Figure 3a,c). Inversely, differences in functional response between the treatment where A. gifuensis was alone and the treatment where A. gifuensis was sharing the experimental patch with H. axyridis were detected in the aggregation or uniform treatment at aphid densities below 150, respectively (Figure 3b,d).  Functional response curves were overlapped between the treatment where H. axyridis was alone and the treatment where H. axyridis was sharing the experimental patch with A. gifuensis for aggregate Insects 2020, 11, 325 7 of 12 and uniform treatments, respectively (Figure 3a,c). Inversely, differences in functional response between the treatment where A. gifuensis was alone and the treatment where A. gifuensis was sharing the experimental patch with H. axyridis were detected in the aggregation or uniform treatment at aphid densities below 150, respectively (Figure 3b,d).

Discussion
In this study, aphid killed rates declined with increasing aphid densities in all cases. Logistic regression analysis indicated that the data from all treatments could fit the type II functional response curve statistically. The type II functional response is common in aphid predators and parasitoids [51,52]. For type II functional response, the unstable enemy-pest dynamic is likely to occur because the predation/parasitism rate would decrease with increasing prey density. Therefore, predators or parasitoids that exhibit type II functional response often cause prey extinction at low densities, but do not affect the prey populations at high prey densities [21,53]. When predators or parasitoids with type II response are applied in biological control systems, a high enemy-pest ratio is necessary to achieve effective pest suppression [2].
In this study, we found that both the killed rates of single H. axyridis and A. gifuensis were affected by prey distribution. Our results are similar to that of Feng et al. [44]. Single H. axyridis or A. gifuensis exhibited higher functional response curves in the aggregate treatment than in the uniform treatment at aphid densities between 40-180 or 70-170, respectively. The possible explanation is that H. axyridis or A. gifuensis may have enough searching time in a 24-h foraging period, enabling H. axyridis or A. gifuensis to encounter, consume or parasitize more aphids when aphid densities were low. With increasing aphid density, both H. axyridis and A. gifuensis may reveal area-restricted foraging behavior in the aggregate treatment like other predators and parasitoids [30,53]. However, when H. axyridis or A. gifuensis are foraging in patches with prey uniformly distributed, they may

Discussion
In this study, aphid killed rates declined with increasing aphid densities in all cases. Logistic regression analysis indicated that the data from all treatments could fit the type II functional response curve statistically. The type II functional response is common in aphid predators and parasitoids [51,52]. For type II functional response, the unstable enemy-pest dynamic is likely to occur because the predation/parasitism rate would decrease with increasing prey density. Therefore, predators or parasitoids that exhibit type II functional response often cause prey extinction at low densities, but do not affect the prey populations at high prey densities [21,53]. When predators or parasitoids with type II response are applied in biological control systems, a high enemy-pest ratio is necessary to achieve effective pest suppression [2].
In this study, we found that both the killed rates of single H. axyridis and A. gifuensis were affected by prey distribution. Our results are similar to that of Feng et al. [44]. Single H. axyridis or A. gifuensis exhibited higher functional response curves in the aggregate treatment than in the uniform treatment at aphid densities between 40-180 or 70-170, respectively. The possible explanation is that H. axyridis or A. gifuensis may have enough searching time in a 24-h foraging period, enabling H. axyridis or A. gifuensis Insects 2020, 11, 325 8 of 12 to encounter, consume or parasitize more aphids when aphid densities were low. With increasing aphid density, both H. axyridis and A. gifuensis may reveal area-restricted foraging behavior in the aggregate treatment like other predators and parasitoids [30,53]. However, when H. axyridis or A. gifuensis are foraging in patches with prey uniformly distributed, they may move more frequently between patches and spend more time in their searching process, and thus decrease aphid predation/parasitism. However, the majority of insect predators are digestion-limited, and the digestion process could affect their foraging efficiency [44]. For H. axyridis, this means that they digest prey slower than they handle them. Therefore, aphid consumption by H. axyridis in two aphid distributions did not differ significantly when aphid densities were above 180. As for A. gifuensis, the parasitism capacity may be limited by the parasitoid egg number in their body like other parasitoids [54], which did not increase parasitism efficiency with aphid density increased to around 170 in either aphid distributions.
Prey distribution could affect the functional response of H. axyridis when A. gifuensis was present. The reasons for the results may be similar to the single H. axyridis treatment. However, the number of aphids parasitized by A. gifuensis sharing the same experimental unit with H. axyridis was not affected by aphid distribution. These results differ from our initial hypothesis. Previous studies found that predators and parasitoids may reveal area-restricted foraging behavior in high prey quality patches [26,33,44], which may increase the antagonistic and intraguild interaction strength between H. axyridis and A. gifuensis. In our experiment, when H. axyridis was present, there was a trend towards increasing the number of aphids parasitized by A. gifuensis in the uniform treatment compared with aggregate treatment, but the differences were not significant. This might be due to the complementary resource use and partition resources by predators and parasitoids [55,56]. In resource partitioning, natural enemy species consume different subpopulations of prey so that a greater proportion of the total prey populations can be exploited by multispecies communities [6,[57][58][59][60]. In the present study, H. axyridis and A. gifuensis maybe partitioning resources across other, unexplored niche axes, and thus did not affect the parasitism efficiency across all aphid densities in the aggregate treatment compared with uniform treatment. Future studies are needed to consider this possibility and explore the mechanisms.
The number of aphids consumed by H. axyridis was not affected by the presence of A. gifuensis compared with single H. axyridis across all aphid densities under two aphid distributions. In fact, the estimated handling times and attack rates of H. axyridis were similar, independently of the presence of A. gifuensis. Previous studies found that female H. axyridis did not exhibit any preference for unparasitized Aphis glycines Matsumura and aphids parasitized by Aphelinus certus Yasnosh [61]. Accordingly, it is possible that H. axyridis did not exhibit a preference between unparasitized M. persicae and parasitized aphids. Nevertheless, parasitoid progeny survivorship was higher in the treatment where A. gifuensis was alone in the aggregate treatment or at aphid densities below 150 in the uniform treatment than in the treatment when A. gifuensis was sharing the same unit with H. axyridis. This means that intraguild interactions occurred between H. axyridis and A. gifuensis. Previous studies found that predators would feed on more prey with increasing levels of prey aggregation because they prefer patches with higher prey densities [30,53]. This may increase the predation rate of H. axyridis on parasitized aphids when all prey were aggregated in one plant. On the aphid uniform treatment, H. axyridis could have enough searching time in the 24-h exposure period when aphid densities were low. With aphid density increasing, H. axyridis may move more frequently between patches and spend more time on foraging prey, and then the encounter probability between the predator and parasitized aphid starts to drop off. This may increase the possibility of A. gifuensis progeny survivorship in the uniform treatment when aphid densities were high.

Conclusions
In the current study, density dependent predation rate of H. axyridis was affected by aphid distribution, but not influenced by heterospecific interactions. The parasitism rate of A. gifuensis was affected by aphid distribution, and by heterospecific interactions in both the aggregate and uniform