The Combined Effects of Warming and Body Size on the Stability of Predator-Prey Interactions

Environmental temperature and body size are two prominent drivers of predation. Despite the ample evidence of their independent effects, the combined impact of temperature and predator-prey body size ratio on the strength and stability of trophic interactions is not fully understood. We experimentally tested how water temperature alters the functional response and population stability of dragonfly nymphs (Cordulegaster boltonii) feeding on freshwater amphipods (Gammarus pulex) across a gradient of their body size ratios. Attack coefficients were highest for small predators feeding on small prey at low temperatures, but shifted toward the largest predators feeding on larger prey in warmer environments. Handling time appeared to decrease with increasing predator and prey body size in the cold environment, but increase at higher temperatures. These findings indicate interactive effects of temperature and body size on functional responses. There was also a negative effect of warming on the stability of predator and prey populations, but this was counteracted by a larger predator-prey body size ratio at higher temperatures. Here, a greater Hill exponent reduced feeding at low prey densities when predators were much larger than their prey, enhancing the persistence of both predator and prey populations in the warmer environment. These experimental findings provide new mechanistic insights into the destabilizing effect of warming on trophic interactions and the key role of predator-prey body size ratios in mitigating these effects.

Environmental temperature and body size are two prominent drivers of predation. Despite the ample evidence of their independent effects, the combined impact of temperature and predator-prey body size ratio on the strength and stability of trophic interactions is not fully understood. We experimentally tested how water temperature alters the functional response and population stability of dragonfly nymphs (Cordulegaster boltonii) feeding on freshwater amphipods (Gammarus pulex) across a gradient of their body size ratios. Attack coefficients were highest for small predators feeding on small prey at low temperatures, but shifted toward the largest predators feeding on larger prey in warmer environments. Handling time appeared to decrease with increasing predator and prey body size in the cold environment, but increase at higher temperatures. These findings indicate interactive effects of temperature and body size on functional responses. There was also a negative effect of warming on the stability of predator and prey populations, but this was counteracted by a larger predator-prey body size ratio at higher temperatures. Here, a greater Hill exponent reduced feeding at low prey densities when predators were much larger than their prey, enhancing the persistence of both predator and prey populations in the warmer environment. These experimental findings provide new mechanistic insights into the destabilizing effect of warming on trophic interactions and the key role of predator-prey body size ratios in mitigating these effects.

INTRODUCTION
The strength of trophic interactions is a key determinant of population, community, and ecosystem stability (McCann et al., 1998;McCann, 2000). Ongoing climate warming may reduce stability, not only by accelerating the risk of species extinctions (Urban, 2015), but also by modulating consumer-resource interactions (Post, 2013;Urban, 2015). Given that environmental warming exacerbates the need for energy uptake by consumers , warmer environments should be characterized by stronger feeding rates of predators up to their thermal tolerance limits. Studies across a large latitudinal gradient support this general prediction, and highlight that arthropods, as ectothermic organisms, are particularly susceptible to the effects of warming (Roslin et al., 2017;Romero et al., 2018). Stronger trophic linkages as a consequence of warming can negatively influence the stability of predator and prey populations (Vasseur and McCann, 2005;Rall et al., 2010;Vucic-Pestic et al., 2011). However, warming can also weaken trophic interactions and put predators at risk of extinction, particularly if their metabolic demands increase more than their predation rates Vucic-Pestic et al., 2011;Fussmann et al., 2014). Coupled changes to interaction strength and energetic efficiency can ultimately destabilize ecological communities at lower temperatures, but have stabilizing effects in warmer environments (Synodinos et al., 2021).
Body size is among the most conspicuous characteristics of all organisms, and it can span up to 20 orders of magnitude from bacteria to whales (Andersen et al., 2016;Blanchard et al., 2017). Body size is often considered to be a "master trait" (Andersen et al., 2016) as it is closely linked to individual metabolism , growth (Gillooly et al., 2002), reproduction (Arendt, 2011), dispersal and locomotion (Jetz et al., 2004). At the community level, body size is a key predictor of population densities and food web structure (Damuth, 1981;Peters, 1983;Cohen et al., 2003;White et al., 2007;Perkins et al., 2019). The relative body size of predators and their prey (typically measured as their body mass) also underlies the strength of trophic interactions and food web stability. In general, experimental studies have indicated a positive association between predatorprey body size ratio and the strength of per capita predation (Emmerson and Raffaelli, 2004). On the other hand, large-bodied predators are disproportionately more prone to extinctions due to human-induced perturbations (Enquist et al., 2020), including overexploitation (Pauly et al., 1998) and climate warming (Petchey et al., 1999). Increasing predation rates in systems characterized by high predator-prey body size ratios and a greater sensitivity of large predators to climate warming can destabilize ecological communities and indicates that both these drivers should be investigated in concert.
Per capita predation rate as a function of prey densitythe functional response of a predator (Holling, 1959(Holling, , 1966)is the main representation of trophic interactions and a central concept in food web ecology. Per capita predation rates commonly increase with prey density in a non-linear manner, with either hyperbolic (Type-II) or sigmoidal (Type-III) relationships (Holling, 1959(Holling, , 1966. These functional response models are characterized by three key parameters. (i) The attack coefficient describes the rate of a successful search through an environmental area or volume. An increase in the attack coefficient translates to steeper ingestion rates at lower prey densities, which can reduce the stability of predator-prey interactions. (ii) Handling time describes the time a predator requires to subdue, ingest, and digest a single prey item. Handling time becomes relatively more important at higher prey densities and increases in this parameter translate to a decrease in maximum ingestion rate. (iii) The Hill exponent (from original work by Barcroft and Hill, 1909) describes how predators  can switch between Type-II and Type-III functional responses.  Smaller values of the Hill exponent tend toward Type-II, whereas  higher values tend toward Type-III functional responses, due  to increased and reduced feeding rates at low prey densities, respectively (Real, 1977).
Several studies have tested the independent impact of warming on functional response parameters. This work indicates that the attack coefficient generally increases, whereas handling time declines with rising temperatures (Rall et al., 2012;Uiterwaal and DeLong, 2018). Fewer experimental studies have investigated temperature effects on the Hill exponent, and those indicated an increase in the Hill exponent (Daugaard et al., 2019;Sohlström et al., 2021) or even a switch between Type-III and Type-II functional responses (Daugaard et al., 2019) after experimental warming. Across a large range of environmental temperatures, there is compelling evidence for a unimodal effect of warming on the functional response parameters, however, pointing to an optimal temperature that maximizes per capita predation rates (Englund et al., 2011;Uszko et al., 2017;Uiterwaal and DeLong, 2020).
The attack coefficient has been suggested to exhibit a humpshaped relationship with predator-prey body mass ratio McCoy et al., 2011;Rall et al., 2011), underpinned by increasing attack coefficients for larger predators and weaker or no effects of prey size (Rall et al., 2012;Uiterwaal and DeLong, 2020). In contrast, handling time has a U-shaped relationship with predator-prey body mass ratio (Rall et al., 2012), underpinned by decreasing handling time for larger predators and smaller prey Kalinkat et al., 2011). Taken together, these findings suggest that there is an optimal predator-prey body size niche for foraging, whereby larger predators generally have a higher chance to encounter, successfully attack, and handle a single smaller prey (Gergs and Ratte, 2009), but the prey should become harder to perceive if they are too small, and harder to handle if they are too large. A recent meta-analysis has also suggested that separate body masses of predators and prey are better at explaining predation rates than the combined predator-prey mass ratio (Uiterwaal and DeLong, 2020). These findings have led to the development of a generalized allometric functional response framework, which allows a continuous transition between Type-II and Type-III functional responses depending on predator and prey body masses, and more accurate modeling of predator-prey dynamics (Kalinkat et al., 2013).
Although independent effects of warming and body size of predators and prey on the strength and stability of trophic interactions have received substantial research attention, much less is known about their combined impacts. This poses a critical gap in our knowledge as global climate change is not only increasing environmental temperatures (IPCC, 2018), but also altering the body size of predators and their prey . The temperature-size rule and temperature-induced changes in community composition led to the suggestion that declining body mass may be a universal response to climate warming (Atkinson, 1994;Daufresne et al., 2009;Evans et al., 2020), although there are some emerging examples of the contrary (Gardner et al., 2011;Ruger and Sommer, 2012;O'Gorman et al., 2017). Smaller organisms can change their body size and shape in response to warming more rapidly than larger organisms due to their higher mass-specific metabolic rates and shorter generation times (Gillooly, 2000;Savage et al., 2004;Tan et al., 2021). This would result in differential temperature-size responses across trophic levels and an increase in predator-prey body mass ratios under the future warmer climate. A recent analysis of marine organisms also showed that warming increases the mean body size of prey consumed by smaller predators but reduces the size of prey consumed by larger predators , with potentially strong and complex consequences for the strength and dynamics of trophic interactions. Thus, it is vital to mechanistically understand both the independent and combined impacts of temperature and body mass on functional responses, given the likelihood that these two factors will play an increasingly important role in altering trophic interactions in our changing world.
Here, we experimentally compared functional responses of a model predator-prey pair at three temperatures and across a gradient of predator-prey body mass ratios to determine the combined effects of warming and body size on the components of predation. Based on two recent syntheses of functional response experiments, we hypothesized (H1) that warming and increasing predator and prey body mass should increase the attack coefficient (Figures 1A,B; Rall et al., 2012;Uiterwaal and DeLong, 2020). We also expected (H2) that warming and increasing predator body mass and decreasing prey body mass should reduce the handling time (Figures 1C,D; Rall et al., 2012;Uiterwaal and DeLong, 2020). Finally, we predicted (H3) a decline in the Hill exponent, or an outright switch from a stabilizing Type-III to destabilizing Type-II functional response (Oaten and Murdoch, 1975), with increasing temperature and prey body mass and decreasing predator body mass (Figures 1E,F; Kalinkat et al., 2011;Daugaard et al., 2019). If the impacts of temperature and body mass are independent, we would only expect warming to alter the intercept or strength of the relationship (i.e., a change in the slope, but with the same sign) between body mass and functional response parameters (Figure 1). Alternatively, if there are interactive effects of temperature and body mass, warming should alter the directionality (i.e., a change in the slope and the sign) or shape of the relationship (i.e., a shift to a hump-shaped relationship). To test whether and how these responses to warming and body size influence persistence and stability of trophic interactions, we used the experimental data to parameterize a bioenergetic predator-prey model and simulated the population dynamics of both trophic levels.

Experimental System: Predator and Prey
To test our hypotheses, we used model predator and prey species that are commonly found in streams throughout Europe.
Cordulegastridae are freshwater dragonfly predators typically associated with lotic waters, streams, and brooks, where their nymphs often inhabit leaf packs and sediment. Cordulegaster boltonii (the golden-ringed dragonfly) is the most abundant Cordulegastridae species in Europe, ranging from southern Spain to southern Scandinavia (Askew, 1988). The species has 14 instars during larval development, which range from 0.4 mm up to more than 8.0 mm of head capsule width (HCW) (Woodward and Hildrew, 2002a). Cordulegaster boltonii is a sit-and-wait predator (Woodward and Hildrew, 2002b) with one of the fastest attack responses among all animals (Pritchard, 1965). It has a wide diet, ranging from small to very large prey which can be captured by the use of their labial mask and palps (Corbet, 1999). The model predators were hand-collected from Broadstone Stream, SE England, where there was an eruption of C. boltonii in the mid-1990s, with densities reaching 72 nymphs m −2 , increasing top-down control of the stream community (Woodward and Hildrew, 2001).
Gammarus pulex is an amphipod species which inhabits a wide variety of habitats due to its high adaptability to light availability, temperature fluctuations, and food sources. Gammarus pulex is an herbivorous shredder and detritivore (Friberg and Jacobsen, 1994;Rong et al., 1995;Gayte and Fontvieille, 1997), but there is also evidence that it feeds on fungi and algae which grow on detritus (Moenickes et al., 2011) and it can occasionally behave as a predator (MacNeil et al., 1997) and a cannibal (Summers et al., 1996;Kelly et al., 2002). The species plays a key role in many stream food webs by shredding leaf material, making it available for decomposers and filter feeders in the stream (Cummins and Klug, 1979), thus linking energetic transport from terrestrial to freshwater environments (Cummins and Klug, 1979;Graca et al., 2001;Felten et al., 2008). The life span of the amphipod varies between 1 and 2 years, with the highest growth rates recorded shortly after birth (Sutcliffe et al., 1981). Gammarus occur in Broadstone Stream in low densities, so the model prey for this study were hand-collected from the River Cray and several streams in the greater London area, where they reach high densities.
Functional response experiments were conducted in May-July of 2013 and 2016 in constant temperature rooms located at Queen Mary University of London, United Kingdom. The attack coefficient is predicted to increase with predator body mass, with a much weaker increase for prey body mass, while the intercepts and slopes of these relationships are expected to increase with warming. (C,D) Handling time is predicted to decrease with predator body mass and increase with prey body mass, while the intercepts and slopes of these relationships are expected to decrease and increase with warming, respectively. (E,F) The Hill exponent is predicted to increase with predator body mass and decrease with prey body mass, while the intercepts and slopes of these relationships are expected to decrease with warming.
The rooms were set to 5, 10, and 15 • C (which we will refer to henceforth as cold, tepid, and warm environments, respectively) and to a 12:12 h light:dark cycle. Wild populations of freshwater macroinvertebrates are known to be adapted to local environmental temperatures (Maazouzi et al., 2011) and the targeted range of experimental temperatures is within the range experienced by amphipods and odonates in United Kingdom rivers (Hannah and Garner, 2015). Identical 2-L plastic aquaria (18 cm long × 11.7 cm wide × 11 cm deep, SAVIC, Belgium), containing 2 cm of washed gravel substrate served as experimental arenas. Each experimental unit consisted of one predator individual and prey of different densities: 1, 2, 3, 4, 8, 16, 32, 48, or 64 individuals per experimental arena. We randomly assigned treatment combinations of predator body size, prey body size, and prey density in each constant temperature room. Predators were starved for at least 2 days prior to the beginning of each experiment. They were placed in the arenas for 1 h prior to the introduction of prey to allow them to acclimatize, which requires time to bury itself in the sediment for camouflage. After 24 h, predators were removed from each arena. The water and gravel were then sieved and sorted by hand to count the remaining prey individuals. Between three and five replicates were conducted for each treatment giving a total of 884 experiments.

Model Fitting and Parameter Estimates
We fitted three different functional response models (Type-II,  strict Type-III, and generalized Type-III) to the experimental data. The Type-II functional response can be described as: where F is the feeding rate, N is the density of the prey, a is the attack rate, and h is the handling time. The Type-III functional responses can be described as: where b is the attack coefficient and 1+q is the Hill exponent. Note that q = 1 for the strict Type-III functional response, whereas q is flexible (though always > 0) for the generalized Type-III functional response. We incorporated body mass into the functional responses using established allometric scaling of the attack rate, attack coefficient, handling time, and Hill exponent as follows: where m c is the consumer (i.e., predator) mass, m r is the resource (i.e., prey) mass, and R = m c m r is the predator-prey body mass ratio. h 0 , a 0 and b 0 are normalization constants, α r and β r are scaling exponents for resource mass, α c is a scaling exponent for consumer mass, ε is a scaling exponent for predator-prey body mass ratio, and q max denotes the maximum of the sigmoidal shaped curve for q, where the half-saturation density q 0 defines the predator-prey body mass ratio for which q = q max × 0.5.
All model fitting and analyses were performed in R 4.0.2. All three functional response models (Type-II, strict Type-III, and generalized Type-III, with a or b, h and q predicted by the allometric model) were fitted to the three temperature treatments separately. We conducted numerical simulations of the following ordinary differential equation (ODE): These numerical simulations were combined with Bayesian parameter estimation by drawing samples from the posterior probability distribution P(θ| y) of the model parameters θ (scaling constants and exponents of Equations 4-7) given the data y, based on the likelihood function P(y| θ) and prior distributions P(θ). We coded the models using the "rstan" package (Stan Development Team, 2020). In each iteration of the Markov Chain Monte Carlo (MCMC) sampling (for a current sample θ), the numerical solution of Equation 8 was computed with the built-in Runge-Kutta ODE solver, to predict prey densitieŝ N end at the end of the experiment for each given initial prey density N start . The likelihood was evaluated assuming binomial distributions of the observed numbers of eaten prey N eaten with n = N start trials and p = N start −N end N start success probabilities. We used normal distributions with zero means, but different standard deviations (SD) for all priors. For the scaling exponents α r , α c , β r , and ε we used weakly informative priors in the range of previous results (SD = 1, Kalinkat et al., 2013), but vague priors (a wider SD = 10) for the normalizing constants a 0 , b 0 , and h 0 . For the generalized Type-III, we used weakly informative priors for the exponent maximum q max (SD = 10) and the halfsaturation density q 0 (SD = 1,000, corresponding to the range of experimental predator-prey body mass ratios).
We fitted all models by running four individual MCMC chains with an adaptation phase of 1,000 iterations and a sampling phase of 2,500 iterations each, summing up to 10,000 samples of the posterior distribution. Visual inspection of the trace plots and density plots showed a good mixture of the chains. Gelman-Rubin statistics ofR < 1.01 and an adequate effective sampling size (i.e., the estimated number of independent samples) verified convergence (Gelman and Hill, 2006). The Leave One Out Information Criterion (LOOIC) was used for model comparison, which was computed from the log-likelihood values of the posterior samples in the "loo" package (Vehtari et al., 2017(Vehtari et al., , 2020.

Stability Analysis
The best performing functional response model was used to parameterize a Rosenzweig-MacArthur population dynamics model (Rosenzweig and MacArthur, 1963;Williams et al., 2007) for three temperatures and three predator × three prey body size combinations: Here, the prey population N [individuals arena −1 day −1 ] undergoes logistic growth with population growth rate r and carrying capacity K, offset by a population-level predation rate, F(N)P, based on the best-fitting functional response model. The predator population P [individuals arena −1 day −1 ] increases based on the same functional response scaled by energetic efficiency ω and predator-prey body mass ratio R, offset by a loss term defined by metabolic rate I. We modeled the temperature dependence of population growth rate for G. pulex as per Sutcliffe et al. (1981): where T is environmental temperature [ • C]. We determined the temperature dependence of carrying capacity (see Supplementary Figure 1) using seasonal changes in abundance of G. pulex in an English stream from Welton (1979): where ω 0 [e 2.266 ] is a normalization constant, E ω [0.164 eV] is an activation energy, and T 0 [293.15 K] was the normalization temperature used in T Arr . Finally, the temperature and body mass scaling of metabolic rate for the predator was determined from laboratory experiments on dragonfly nymphs (see Supplementary Figure 2):  (Cummins and Wuycheck, 1971). For each simulation, we defined stability as: where the range of predator densities was measured after the system reached a steady state or limit cycle. In this way, stability = 1 for steady states and declines toward 0 with an increasing amplitude of predator densities.

RESULTS
The most parsimonious model was the generalized Type-III functional response (LOOIC = 3,658.8, r 2 = 0.63) across all treatment combinations, which had a considerably lower LOOIC value than both the Type-II model ( LOOIC ± SE = 69.6 ± 25.8, r 2 = 0.62) and the strict Type-III model ( LOOIC ± SE = 70.4 ± 37.2, r 2 = 0.62). Parameter value estimates for the best fitting model are shown in Table 1, with model fit to the data illustrated in Figure 2.

Attack Coefficient
The highest attack coefficients were found for the smallest predators feeding on the smallest prey in the cold environment (Figures 3A,D,G), but this shifted toward the largest predators consuming medium and large prey in the tepid (Figures 3B,E,H) and warm (Figures 3C,F,I) environments. This offers partial support for our first hypothesis, with attack coefficients generally increasing with predator body mass as temperature also increased, though with opposing effects for prey body mass. Importantly, this also shows an interactive effect of temperature and body size, with temperature altering the shape and directionality of the relationship between predator and prey body mass and the attack coefficient (compare with Figures 1A,B)

Handling Time
The lowest handling times occurred in the cold environment, where they also tended to decrease with increasing predator and prey body mass (Figures 4A,D,G). In contrast, handling times were generally greater and increased with both predator and prey body mass in the tepid (Figures 4B,E,H) and warm (Figures 4C,F,I) environments. This offers partial support for our second hypothesis, with handling times declining with decreasing prey mass as temperature also increased, though with opposing effects for predator body mass. It also indicates an interactive effect of temperature and body size, with temperature potentially altering the directionality of the relationship between predatorprey body mass and handling time (see Figures 1C,D). Parameters b 0 , β r , and ε are the normalization constant and body mass scaling exponents for the attack coefficient; h 0, α r , and α c are the normalization constant and body mass scaling exponents for handling time; while q 0 and q max define the shape of the curve produced by the Hill exponent.
FIGURE 2 | Generalized Type-III functional response curves fitted to the experimental data for each temperature × predator body mass × prey body mass combination. Panels are arranged from (A-I) in approximately increasing order of predator-prey mass ratio. Solid blue, dashed yellow, and dotted red curves (with transparent areas of the same color) represent the model fit ± 95% CI in the cold (5 • C), tepid (10 • C), and warm (15 • C) temperature treatments, respectively. Solid symbols (blue squares, yellow circles, and red triangles) represent individual feeding experiments.

Hill Exponent
The largest Hill exponents occurred in the cold environment and appeared to decrease with warming ( Figures 5A-C), although the effect of temperature was associated with large 95% posterior probability intervals (Figure 5). The Hill exponent also decreased with increasing prey body mass (Figures 5D-F) and decreasing predator body mass (Figures 5G-I). These findings offer support for our third hypothesis, though note that there was no evidence for an outright switch from Type-III to Type-II functional response with warming or changing predator and prey body mass ( Table 1). These results also indicate independent effects of temperature and body size on the Hill exponent, with temperature only altering the strength and not the directionality or shape of the relationships (see Figures 1E,F).

Population Stability
Predator-prey population dynamics were predicted to be stable for all body mass ratios in the cold environment  Figure 4), but became much more variable in both the tepid (Supplementary Figure 5) and warm (Supplementary Figure 6) environments. An increase in predator-prey body mass ratio (i.e., when predators were much bigger than their prey) led to less variable population dynamics, however, and thus dampened the destabilizing effect of temperature (Figure 6). Furthermore, the most variable dynamics and thus lowest stability occurred at intermediate body mass ratios (Figure 6), indicating a non-linear effect of predator-prey body mass ratio on population stability in warmer environments.

DISCUSSION
Our findings demonstrate interactive effects of warming and body mass on the mechanistic components of predation, particularly attack coefficients and handling times. Although previous studies have identified many similar independent effects of body size and temperature across the experimental gradients we used here (Englund et al., 2011;Rall et al., 2012;Uiterwaal and DeLong, 2020), we show that these factors need to be considered in tandem (i.e., interactively) when interpreting functional response data or modeling trophic relationships in the face of ongoing global warming. The Hill exponent also declined with prey body mass and increased with predator body mass, with no effect of warming on these relationships. Although not universal, larger Hill exponents are often associated with increased stability due to reduced predation at low prey densities (Uszko et al., 2015;Daugaard et al., 2019). This meant that the destabilizing impacts of warming on population biomasses were alleviated by an increase in the predator-prey body mass ratio and the associated larger Hill exponents.
In partial agreement with our first hypothesis, the attack coefficient tended to increase with predator body size, but only when they fed on medium and large prey in the warmer environments (Figures 3H,I). Large predators were less efficient at successfully attacking small prey individuals, which is in line with findings across a wide range of taxa showing that large and small consumers preferentially feed on large and small resources, respectively (Jonsson et al., 2005;Brose et al., 2006). This was particularly clear in the cold environment, where the largest attack coefficients occurred when small predators fed on small prey. Large predators exhibited a very low rate of successful attacks in the cold environment, indicating that they may have been largely inactive and require higher temperatures to start foraging. It also suggests that these predators are likely to have very low foraging rates during the colder winter months, with their strongest interactions occurring in summer, generating strong seasonal differences in energy flow through the food web.
Our findings also provide strong support for a hump-shaped relationship between attack coefficient and predator-prey body mass ratio, which has been widely reported in the literature McCoy et al., 2011;Rall et al., 2011). This indicates that C. boltonii has an optimum prey size that can be successfully attacked, and they may not perceive prey that are > 1,000 times smaller, or struggle to subdue prey that are < 10 times smaller than themselves (Supplementary Figures 3A-C). Environmental temperature altered the directionality of the relationship between handling time and body size of predators and prey, which has not previously been demonstrated in functional response experiments. This interactive effect between temperature and body size should be treated with caution because highest prey density treatments were logistically difficult to establish, which resulted in either low or no replication (see Supplementary Table 1). However, all nine mass ratio treatment combinations were fitted together, with these handling times being informed by the data across all mass ratio combinations (Equation 6), increasing the robustness of these parameter estimates. Overall, the effects of prey body size were stronger (i.e., steeper slopes) than the effects of predator size. The predicted positive relationship between handling time and prey body mass (Rall et al., 2012;Uiterwaal and DeLong, 2020) was reversed in the cold environment, where predators were more efficient at handling large prey. These large and energetically important prey individuals are likely to evade the predators in tepid and warm environments, but their slower movements might make them easier to catch and subdue in cold environments. Interestingly, the predicted negative relationship between handling time and predator body mass (Rall et al., 2012;Uiterwaal and DeLong, 2020) was only observed in the cold environment. Although our experimental temperatures were well within the normal range over which such relationships are investigated (Englund et al., 2011;Uiterwaal and DeLong, 2020), longer handling times in warmer environments may indicate that predator foraging performance is impeded by warming, or that their prey are more actively trying to escape them (Domenici et al., 2019).
The Hill exponent appeared to decline with experimental warming, which adds to the broad range of contrasting responses that have been described recently (Uszko et al., 2017;Daugaard et al., 2019;Sohlström et al., 2021). However, this effect was associated with wide uncertainty and should be interpreted with caution. The associated high feeding rates at low prey densities most likely underpinned the negative effect of warming on the stability of predator and prey populations (Figure 6). Similar negative warming-stability relationships have recently been described at low and moderate environmental temperatures (Synodinos et al., 2021). Interestingly, our bioenergetic model simulations indicate that this destabilizing effect was buffered by larger predator-prey body mass ratios through the independent effects of predator and prey body mass on the Hill exponent. We found the greatest Hill exponents when the largest predators were feeding on the smallest prey individuals, which is similar to the stabilizing switch from Type-II to Type III functional response observed at the highest predator-prey body mass ratios in previous research . In combination with the high handling time for large prey in the warm environment, this may indicate a greater feeding efficiency for predators when consuming highly abundant small prey that are easier to catch and subdue. Thus, warming may alter the optimal size niche of predator-prey interactions by shifting the risk-reward tradeoff between energy expenditure (i.e., foraging) and energy intake (i.e., successful predation). Future evidence for such a trade-off could come from an experimental characterization of the thermal performance curves for predator and prey species, coupled with mechanistic modeling of their population dynamics.
The lowest stability of predator and prey populations occurred at intermediate body mass ratios, suggesting a non-linear relationship between the two variables. While it may seem counterintuitive that the highest body mass ratios led to the greatest stability (e.g., in contrast to Jennings and Warr, 2003), our manipulated predator-prey size ratios only spanned the range of values where trophic interactions were still feasible (i.e., the largest and smallest predators could still consume the smallest and largest prey, respectively). Expanding this range to include much smaller or larger prey should eventually destabilize the trophic interactions through the inability of the predator to consume the prey. Overall, the application of a generalized allometric functional response (Kalinkat et al., 2013) allowed us to discover this novel stabilizing mechanism of larger body mass ratios under environmental warming, and it should be utilized more widely in predator-prey interaction research.
It is important to note that our experiments were conducted with just one model predator-prey pair. Further research is thus needed to determine the generality of the observed interactive effects of temperature and body size on functional response parameters and population stability. A growing number of studies also indicate that predator-prey interactions may exhibit different responses depending on the dimensionality of the interacting pair (Pawar et al., 2012;Uiterwaal and DeLong, 2020). Our study involved a sit-and-wait predator feeding on an active prey that can occupy three-dimensional space in the water column. The predator is thus likely to be less affected by warming than active hunters who rely on speed and agility for successful predation events, which are strongly influenced by environmental temperature (Öhlund et al., 2015;Cloyed et al., 2019).
Similarly, interactions involving sessile prey may be more strongly determined by the response of the predator. Particular attention should be paid to investigating the combined effects of warming and body size on functional response parameters across a range of interaction types in future research.
Although previous studies have explored the independent effects of temperature or body size on trophic interactions and stability (Englund et al., 2011;Fussmann et al., 2014;Uiterwaal and DeLong, 2020;Synodinos et al., 2021), our understanding of the interactive effects between these two key variables remains incomplete. Here, we demonstrate how temperature alters the impact of body size on components of predation and how predator-prey body size ratio can buffer the destabilizing effect of warming. We also identify novel mechanisms by which these drivers can alter ecological communities, which should improve predictive models and thus our ability to manage natural systems in the face of climate change. The functional response, formalized by Holling (1959) and recently advanced by renewed interest (Kalinkat et al., 2013;Uszko et al., 2017;Rosenbaum and Rall, 2018), remains a central concept facilitating the mechanistic understanding of food web structure and stability. A growing number of studies now highlight the flexibility of the Hill exponent and its importance for the stability of populations exposed to warming (Daugaard et al., 2019;Sohlström et al., 2021). Further elucidation of how this component of predator-prey interactions responds to environmental stressors such as warming and biological traits such as body size is a promising avenue for future research.

AUTHOR CONTRIBUTIONS
EO'G and PK designed the study. BG and ELH conducted the experiments. BR and EO'G analyzed the data and performed the model simulations. PK wrote the initial draft. All authors edited the manuscript.

FUNDING
This project was funded by internal support from the Queen Mary University of London to PK, BG, and ELH. PK was also supported by the Royal Society, Newton Advanced Fellowship (Grant No. NAF/R2/180791). EO'G was supported by NERC (NE/L011840/1 and NE/I009280/2). BR acknowledges the support of iDiv funded by the German Research Foundation (DFG-FZT 118, 202548816).

ACKNOWLEDGMENTS
We thank Paul Fletcher for his assistance with animal collection and experimental work. We also thank the editor and two reviewers for their insightful comments on an earlier version of this manuscript.