Functional response of Harmonia axyridis preying on Acyrthosiphon pisum nymphs: the effect of temperature

In the current study, we investigated the functional response of Harmonia axyridis adults and larvae foraging on Acyrthosiphon pisum nymphs at temperatures between 15 and 35 °C. Logistic regression and Roger’s random predator models were employed to determine the type and parameters of the functional response. Harmonia axyridis larvae and adults exhibited Type II functional responses to A. pisum, and warming increased both the predation activity and host aphid control mortality. Female and 4th instar H. axyridis consumed the most aphids. For fourth instar larvae and female H. axyridis adults, the successful attack rates were 0.23 ± 0.014 h−1 and 0.25 ± 0.015 h−1; the handling times were 0.13 ± 0.005 h and 0.16 ± 0.004 h; and the estimated maximum predation rates were 181.28 ± 14.54 and 153.85 ± 4.06, respectively. These findings accentuate the high performance of 4th instar and female H. axyridis and the role of temperature in their efficiency. Further, we discussed such temperature-driven shifts in predation and prey mortality concerning prey-predator foraging interactions towards biological control.

The use of integrated pest management strategies based on biocontrol agents has received increasing prominence worldwide 1 . It has been implemented with tremendous success in both fields and greenhouses 2 , in particular with the intent of reducing the large-scale use of pesticides 3 . While the adoption of biological control is desirable, the successful implementation depends upon the comprehensive understanding of predator-prey interactions, owing to their fundamental role towards ecosystem functionality and food web stability 4 . Several methods can be applied to quantifying these interactions 5,6 , including functional response 7 , numerical response 8 , kill rate 9 , and consumption rate 10 . Functional responses describe how the predation rate changes with resource density 11 . Three types of functional responses are typically expected, with consumption rate being linear up to a constant plateau (Type I), parabolic (Type II), or sigmoid (Type III) 12 , depending on whether the parameters of functional response, i.e., the enemy attack rate and prey handling time, vary with prey density. In the type II functional response, the attack rate defines the steepness of the increase in predation with the increase of prey density, and handling time sets the satiation threshold 13 . Natural enemies with high attack rate and low handling time are thought to be the most efficient biocontrol agents 13 . Many sources are known to regulate these functional response parameters 14 .
Temperature is a chief driver of biological systems through the temperature-dependent nature of biological rates (e.g., metabolic rates) 15 . The effect of temperature on biological rates is likely to be realized from physiology up to species level, influencing population growth rates and carrying capacities 16 as well as ecosystem functions. Temperature can affect prey-predator foraging interactions by altering their behavioural or physiological responses 17 . Warming is shown to increase the predator rate of prey consumption by boosting predator's metabolic rates 18 . In order to consume a large amount of prey, a predator should be adept at searching for and handling its prey, so that it may spend more time searching on consumption events than on prey handling attempts, therefore a change of predatory behaviour or functional response may be expected under warming 19,20 . For instance, the handling time decreased, consumption rate increased, and the type of functional response changed from
Logistic regression between the initial aphid densities offered and the proportion of aphid consumed (Na/No) showed all significantly negative values of the linear coefficients P 1 ; indicating a Type II functional response across all growth stages and temperatures tested ( Table 1). The declining consumption with increasing aphid densities (Fig. 5) also confirmed a Type II functional response. The monotonically declining proportion of consumption with increased aphid densities, led to further confirmation of Type II functional responses (Fig. 6). A more linear decline for the proportion of prey eaten with increasing aphid density was noted for 4th instar (Fig. 6d) and female H. axyridis (Fig. 6f) (Fig. 7a), the shortest handling time (T h ) (Fig. 7b), and the maximum predation rate (T/Th) (Fig. 7c) typically at higher temperatures and later growth stages. The attack rates (a) and maximum predation rates (T/Th) were generally low at low temperatures (15 and 20 °C) that started to increase with warming, reaching a peak at 30 and 35 °C. On the other hand, the handling time (T h ) was the lowest at 35 °C and increased with lower temperatures. The estimate of maximum predation rate (T/Th) was the highest for 4th instar H. axyridis followed by the female and male H. axyridis, respectively. The 4th instar and female H. axyridis similarly showed the lower handling times (T h ) than male H. axyridis.

Discussion
Insect ectotherms are known to respond to thermal conditions for their development, biology, and population dynamics 46,47 , as well as through their trophic interactions 48,49 . The current research explored the functional response of H. axyridis foraging on A. pisum at different growth stages, temperatures and evaluated the effect of temperature on host aphid mortality. The temperature ranges we tested are relevant across a range of temperate or sub-tropical regions. Our results showed increasing host aphid mortality with warmer temperatures, meaning warming could lead to faster prey depletion. A low-temperature threshold (i.e., 20 °C) is reported the best  Figure 1. Percent control mortality of Acyrthosiphon pisum nymphs exposed to different thermal conditions under different prey densities. Data are expressed with box plots for temperature (a), prey density effects (b), and through line graphs for density-wise temperature effects (c) on mean percent mortalities (with 95% confidence intervals (CIs)) following significant temperature by density interaction. Box plots are showing the range of data (lower and upper quartiles and extreme values), median, and mean (symbols), and different alphabets indicate significant differences between group means (P < 0.05; Tukey's HSD test). In panel c, the temperature effects are tested within each prey density offered through running independent ANOVAs and significant means are compared according to non-overlapping 95% CI of difference. Symbols '*' , '**' , and '***' denote significance at 0.05, 0.01, and 0.001 levels, respectively, whereas 'ns' indicates non-significant (P > 0.05) differences for ANOVAs performed.   41 , and for other coccinellids, including A. bipunctata, Hippodamia convergens Guérin-Méneville and Coccinella septempunctata L., preying on M. persicae 19,52 , suggesting some generality of this outcome. In our findings, despite its positive effect on aphid consumption, changing temperature did not change H. axyridis functional response type (i.e., from Type II to Type III). A changing functional response type with thermal changes has been reported for many predators 22 , but rarely for coccinellids. A Type II response describes an increase in predation with increasing prey abundance, gradually decelerating to an asymptotic foraging rate at higher abundance 7 . This also generates negative density-dependent mortality, often a characteristic of predators that provide efficient control at smaller resource 53 . Among the three types of functional responses described by Holling 7 , only type producing a positive density dependent mortality is Type III 11 . The only possibility for exhibiting a Type III response in current study is the concentration of predator hunting in high density patches 54 . This mechanism could have operated in the current experiments, however, no evidence of Type III response was found in our data. The other mechanisms for a Type III response include switching behaviour and predator learning that could not have operated in our system, as the experiments were short-term and single prey-based. Furthermore, H axyridis often shows a Type II functional response. Type II functional responses have been found for H. axyridis when preying on immature S. litura 41 , or when preying on C. juglandicola or P. juglandis at different temperatures 40 , Rhopalosiphum padi L. and Sitobion avenae F. with different fertilizer treatments 55 , or different prey types such as Lipaphis erysimi K. (Hemiptera: Aphididae), Cacopsylla chinensis (Hemiptera: Psyllidae) and Danaus plexippus L. (Lepidoptera: Nymphlidae) 42,56,57 , or different growth stages of A. gossypii, M. persicae, Myzus nicotianae S., Aphis glycines Matsumura (Hemiptera: Aphididae) and Diaphorina citri (Hemiptera: Psyllidae) [58][59][60][61] . However, the change of functional response type with respect to prey distribution 62 or prey quality (either pesticide treated or untreated) have been shown for H. axyridis 63,64 , suggesting complex nature of predator-prey interactions, and, emphasizing the need for further assessments of the functional response www.nature.com/scientificreports/ with regard to factors like pesticides 33 . Insecticides from synthetics and biopesticides groups are commonly applied worldwide [65][66][67][68] and have been shown to have profound effects, both positive or negative, on behavioural or physiological responses of predators 39,69 . The attack rate and handling time describe the overall functional response magnitude. The attack rate (also called the space clearance rate or attack efficiency) describes the ability a predator possesses to catch its prey in a given time frame, and handling time describes the time lost from searching per host consumed 70 . A high attack rate means that the predator is adept at quickly removing hosts from the volumes or areas it is searching, and low handling time means how quickly a predator traps, hunts, and digests prey 70 . In our findings, the parameter estimates showed greater variation across predator growth stages, with frequently higher estimates at later growth stages and higher temperatures. Maximum daily predation rates were temperature-dependent, especially their increase, and the handling time of H. axyridis showed an exponential decrease for all growth stages at the lowest thermal conditions, meaning predators will spend less investment on foraging, possibly using that time for resting and decrease predation 71,72 . Conversely, the greater number of prey consumed due to warming may have resulted from a decrease in handling time with a resultant increase in encounter rates 73 . The expected increase in metabolic rates with warming are associated with greater energy demands, which should cause predators to increase food intake and foraging activity. Determining these parameters with respect to phenological stage confirmed poor response by the first three instars as reported earlier 59,74 , whereas increased foraging performance by final instar and adult H. axyridis (especially female when compared with male) 24,45,52 suggesting H. axyridis may provide better biocontrol efficiency at later stages, plausibly owing to better searching efficiencies. The 4th instar requires large meals to attain the required weight for pupation 75 and adult predators have to prepare for reproduction 56 and other functions related to egg maturation or fertilization 76 . Our results of greater predation by the 4th instar and female H. axyridis are supported by previous reports investigating the functional response of coccinellids, namely Adalia tetraspilota (Hope), Hippodamia variegata (Goeze), Harmonia dimidiata (Fab.), and Scymnus syriacus Marseul preying on many aphid species 44,45,77 .
Our results showed strong potential of H. axyridis for A. pisum control. The 4th instar and female H. axyridis emerged as the best performing biocontrol candidates. We determined the strong role of temperature in predator's efficiency, as accelerated predator action under warming (30 and 35 °C) increased prey consumption. However, there was no evidence for functional response type change with warmer temperatures. Another striking result was that prey mortality increased with warming. This indicates that warming may reduce prey availability by increasing predator action and enhancing prey mortality, which allows us to suggest that H. axyridis release at warmer temperatures may not be feasible against this aphid pest. The 4th instar and female H. axyridis can be    Table 1. Results of the logistic regression analysis of the proportion of nymphs of Acyrthosiphon pisum predated by all stages of Harmonia axyridis relative to the initial number of nymphs provided. www.nature.com/scientificreports/ best used at low temperatures (between 20 and 25 °C) against this aphid pest. Prey depletion under warming can affect predator-prey interactions by modifying functional responses that can then destabilize communities through intraguild predation or other antagonistic interactions triggered in response to prey depletion 28,78 . This necessitates the need for further studies exploring alternative prey resources to support this generalist predator in the management of aphid pests under warmer conditions.

Methods
Cucumber (Cucumis sativa L., cv. Negin; Cucurbitales: Cucurbitaceae) and broad bean (Vicia faba L.; Fabales: Fabaceae) seedlings were grown from seeds purchased from Caoxian County, Shandong, China. The seeds were sown in pots (25 cm diameter, 30 cm deep) with (3:1) soil: manure. The seedlings were maintained under greenhouse conditions of 12-26 °C, 45-55% RH, and 16:8 h (Light: Dark) photoperiod, and subsequently used for rearing and conducting functional response assays. The plant materials used were obtained with prior permission, and the present study is in compliance with relevant guidelines and legislation. For establishing A. pisum culture, the initial populations of aphid collected from unsprayed alfalfa fields were subsequently brought to the laboratory and reared on broad bean plants inside net cages (20 × 10 × 30 cm height). The stock culture of H. axyridis was developed from a pre-established laboratory colony, already available in the same laboratory. The predator was reared on A. pisum infested bean plants (7-8 leaves) inside net cages (60 × 42 × 30 cm height) for three consecutive generations at laboratory conditions of 24 ± 1 °C, 65 ± 5% RH and 16:8 h (Light: Dark) photoperiod. Bean plants were checked daily for predator eggs. When found, egg batches were carefully removed, placed on tissue paper in Petri dishes (9 cm), and transferred to a computer-operated growth chamber, maintained at settings of 25 ± 1 °C, 65 ± 5% RH and 16:8 h (Light: Dark) photoperiod. The www.nature.com/scientificreports/ post-emergence larvae were separated and reared in Petri dishes containing aphids as their diet, refreshed daily. The whole culture was maintained at the Department of Plant Protection, Huazhong Agricultural University, China.
The experimental arena consisted of clear Petri dishes (9 cm diameter), with a micromesh screen over the top for ventilation and bottom covered with clean cucumber leaf disk. The desiccation of cucumber leaf disc was prevented by adding 1% agar solution 79 . The assays were performed with H. axyridis larvae (i.e., 1st instar, 2nd instar, 3rd instar, 4th instar) and adults (male, female) at five constant temperatures (i.e., 15, 20, 25, 30, 35 °C). The homogeneity of predator age was maintained within each tested growth stage. The first instar larvae were separated one by one shortly after hatching to avoid sibling cannibalism. Hatchlings were reared in Petri dishes (9 cm diameter) until maturity on 4th instar nymphs (100-150 aphids/day). Female H. axyridis included mated individuals 59 . First instar larvae were starved for about 6 h, whereas subsequent instars/stages were starved for 24 h to standardize hunger level, according to Islam, et al. 41 . The moist cotton roll offered humidity to all predators during starvation. The use of 4th instar aphid was ensured throughout the experiments as a way to prevent predator preference switch according to prey size 80 . Using a fine camel hairbrush, the aphids at different densities (i.e., 2,4,8,16,32,64,128, and 160 aphids) were transferred in Petri dishes, allowed to spread and settle over the substrate for 30 min, and thereafter were transferred to a computerized growth chamber at the experimental temperatures (i.e., 15, 20, 25, 30, 35 °C), with 70 ± 5% RH and a 16:8 h (Light: Dark) photoperiod. The whole experiment was replicated 10 times for each prey density, growth stage, and temperature combination. The numbers of aphid consumed were recorded every 24 h. Five control replicates were performed for assessing the prey mortality concerning thermal conditions imposed and prey densities offered. Control replicates were kept free from H. axyridis to account for natural mortality and to correct for A. pisum consumption by the predator as a function of natural mortality. Predation mortality data were corrected for control mortality by applying Abbott's correction 81 .
The control mortality data were analyzed between temperatures, aphid densities, and their interaction, by using Univariate Analysis of Variance (ANOVA) in SPSS (version 21), fitting the above three variables as fixed factors against the dependent variable (i.e., host mortality). Significant (P < 0.05) effects were further compared by using Tukey's Honestly Significant Difference (HSD) multiple comparisons test. Prior to analysis, the mortality data were tested for normality and homogeneity of error variance (i.e., homoscedasticity) using Shapiro-Wilk and Levene tests, and Y = √ x + 1 transformed to improve compliance with these assumptions. All means and standard errors in text and figures are calculated with untransformed data.
Aphid consumption by H. axyridis for temperature, growth stage, density, and their two-way and threeway interactions were analyzed by using Generalized Linear Models (GLM) in SPSS (version 21). Kolmogorov-Smirnov test confirmed non-normal distributions of data (P > 0.05), and due to over-dispersion, the data were fitted with negative binomial distribution and a log link function, and factors and interaction effects were analyzed by using the Wald Chi-Square test for a confidence level (CI) of 95%. If needed, the multiple follow up tests (GLM) were run to analyze the temperature and growth stage effects separately across prey densities offered, using SPSS (version 21), and the significance for each test was adjusted by following Bonferroni correction (i.e., dividing the standard P-value criterion by the number of tests) to avoid Type 1 error.
Analysis of the functional response was done in two different phases 11 , in the R statistical environment 82 . The first phase involved the determination of type and estimation of the parameters of the functional response curve. It is compulsory to find the type of functional response for calculating the functional response parameters using a proper model. The type was determined by applying logistic regression of the proportion of prey eaten as a function of initial prey density offered. A polynomial logistic regression equation assuming a binomial distribution of data to define the type of functional response 11 (Eq. 1) was fitted as under: where N a and N o indicate the number of prey consumed and the initial prey density offered, respectively, and Na No is the proportion of prey consumed. The P o , P 1 , P 2 , and P 3 are the regression parameters representing intercept or constant, linear, quadratic, and cubic coefficients, respectively. The coefficients were calculated using maximum likelihood. The values of the linear and quadratic coefficients indicate the nature of functional response, either Type II or Type III. When the value of a linear parameter is negative, the functional response is Type II, and if it is positive with a negative quadratic coefficient, then response is of Type III. The Type II response shows that the proportion of prey consumption decreases as the prey density increases, and a Type III response represents that the proportion of prey consumed increases until an inflection point and then decreases 11 . Once the functional response type was determined, the second phase started where functional response parameters were determined. Data were fitted to the Rogers' type II random predator equation, using non-linear least square regression, as the prey was not replaced during the entire experiment 83 . The attack rate (a) and handling time (Th) were calculated by using the random predator model as under (Eq. 2): where N a is the number of prey eaten, N o is the initial prey density (prey.arena −1 ), a is the attack rate (arena. hour −1 ), T h is the handling time (hour.prey −1 ), and T is time available for predator during the experiment (here 24 h). Here, the "glm" function was used to fit the logistic regression, and the parameters (attack rate a and handling time T h ) of functional response were estimated by using FRAIR (Functional Response Analysis in R, version 4.0.0) 84 . The maximum theoretical predation rate per day (K = T/T h ) 85 indicates the maximum amount of prey that a predator can consume in a given time frame (here 24 h).