Evolution of predator-prey interactions during range expansion in an aquatic insect predator

3060 | Functional Ecology. 2022;36:3060–3072.wileyonlinelibrary.com/journal/fec Received: 4 April 2022 | Accepted: 12 October 2022 DOI: 10.1111/1365-2435.14208 RESEARCH ARTICLE Evolution of predator– prey interactions during range expansion in an aquatic insect predator José A. Carbonell1,2 | YingJie Wang1 | Arnaud Sentis3 | Robby Stoks1 This is an open access article under the terms of the Creative Commons AttributionNonCommercial License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited and is not used for commercial purposes. © 2022 The Authors. Functional Ecology published by John Wiley & Sons Ltd on behalf of British Ecological Society. 1Evolutionary Stress Ecology and Ecotoxicology, University of Leuven, Leuven, Belgium 2Department of Zoology, Faculty of Biology, University of Seville, Seville, Spain 3INRAE, AixMarseille Université, UMR RECOVER, AixenProvence, France Correspondence José A. Carbonell Email: [email protected] Funding information Fonds Wetenschappelijk Onderzoek, Grant/Award Number: G.0524.17N and G.0956.19N; KU Leuven, Grant/Award Number: C16/17/002; European Union; Ministry of Education of Taiwan Handling Editor: Caroline Williams Abstract 1. Many ectotherms are shifting their distributions polewards, which has been associated with the evolution of phenotypic traits and their thermal plasticity. Trophic interactions may determine the dynamics and ecological impact of range expansions. However, it is largely unknown how trait evolution in edge populations shapes trophic interactions. 2. We studied evolutionary changes in the shortterm (functional response) and longterm predator– prey interactions between an aquatic insect predator (the damselfly Ischnura elegans) and its prey (the water flea Daphnia magna) during the predator's ongoing poleward range expansion in northern Europe. 3. Using a commongarden warming experiment at 20 and 24°C we tested for differentiation between predator populations from edge and core regions in metabolic rate and functional response parameters, and used these empirical data to estimate the effects of range expansion on the shortand longterm predator– prey interaction strengths. 4. Metabolic rates did not differ between populations from edge and core regions nor between rearing temperatures. Functional response parameters and their thermal plasticity showed signals of evolution during the range expansion. Attack rates did not differ between predators from edge and core regions, but only decreased under warming in predators from the edge region. Handling times decreased under warming in predators from the edge region but increased under warming in predators from the core region. While handling times were shorter in predators from the core region at 20°C, these did not differ between regions at 24°C. As a result, the shortterm interaction strength was higher for predators from the core region at 20°C, but not different between regions at 24°C. The predator– prey system from the edge region showed lower longterm system stability at 20°C, but this region difference disappeared under warming because the edge region stability then increased. 5. Our results suggest that rapid evolution of functional response parameters during a predator's range expansion reduced the direct feeding impact on its prey and made the predator– prey system from the edge region more unstable, | 3061 Functional Ecology CARBONELL et al. 1 | INTRODUCTION Many ectotherm species are shifting their distributions in a poleward direction in response to climate change and thereby may encounter novel environmental conditions to which they need to adapt (Spence & Tingley, 2020), including novel thermal regimes (Lancaster et al., 2015). It is therefore crucial to study the evolution of species traits associated with range shifts and the evolution of thermal plasticity of those traits to mechanistically understand current range dynamics and inform future ones (Diamond, 2018; Nadeau & Urban, 2019). Accordingly, an increasing number of studies have documented rapid evolution and thermal plasticity of life history traits during range expansions (e.g. Carbonell et al., 2021; Diamond, 2018; Stuart et al., 2019). Yet, range expansion dynamics and their ecological impact may critically depend on interactions with other species and the effects of thermal changes on these interactions (Aguilera et al., 2019; Jones et al., 2020; Nadeau & Urban, 2019). Nevertheless, how trait evolution during range expansions modulates trophic interactions has been rarely studied. Trophic interactions, and particularly predator– prey interactions, can strongly affect the structure, stability and dynamics of aquatic food webs (Batzer & Sharitz, 2006; Paine, 1966), yet how the trophic interaction strength between predators and prey changes during range expansions is largely unknown. Although few studies demonstrated longterm evolution of predator– prey trophic interactions between latitudes (Barton, 2011; De Block et al., 2013; Tran et al., 2016; Wang et al., 2021), addressing changes in trophic interactions during ongoing range expansions would provide key insight into their contemporary, rapid evolution. One notable exception comparing core and edge populations is a mesocosm study showing that larvae of the dragonfly Crocothemis erythraea from an edge population had a lower impact on zooplankton densities compared to those from three core populations, but only at a higher temperature (Therry et al., 2019). While this study nicely illustrated the importance of considering thermal plasticity when studying the ecological effects of range expansions, it did not explicitly study predator traits directly linked to the predator– prey interaction strength. In the shortterm, the interaction strength can be directly estimated based on the functional response, which describes the relationship between per capita predation rates and prey densities. The shape and type (I, II or III) of the functional response curve is determined by two parameters, the predator attack rate that determines feeding rate at low prey densities, and the predator handling time that determines the maximum predation rate at high prey densities (Holling, 1959). The ratio of both parameters is a measure of the shortterm interaction strength (Cuthbert et al., 2019). Under warming, attack rates typically increase, while handling times shorten (Bideault et al., 2021; Thompson, 1978; Twardochleb et al., 2020; Wang et al., 2021), leading to a higher shortterm interaction strength. Yet, thermal responses may be more complex and take the form of unimodal relationships with attack rates decreasing again at very high temperatures when species are near their upper thermal limit (Englund et al., 2011; Sentis et al., 2012). Despite the importance of the functional response to understand trophic interactions, its rapid evolution has been rarely studied in natural systems (but see Urban et al., 2020), and never in the context of range expansions. Given the predicted evolutionary changes in energetically costly life history and dispersal traits during range expansion (Fitt et al., 2019) ask for a greater prey consumption, this is expected to result in the evolution of higher attack rates and/or shorter handling times (Urban et al., 2020). While the functional response may give important insights in the shortterm stability of predator– prey interactions (Hassell et al., 1977), in the longterm also numerical responses across generations of the predator and the prey have to be taken into account (Fussmann et al., 2014). Integrating the functional response parameters with, among others, predator metabolic rates and prey demographic parameters allows us to estimate the longterm interaction strength, the effect of the predator population on the prey population density (Fussmann et al., 2014; Rall et al., 2010). The longterm interaction strength is expected to decrease under warming. This is mainly because warming may cause higher prey intrinsic growth rates that increase the predator equilibrium densities, and higher predator metabolic rates that limits predator population density and thus increase the prey equilibrium densities (Fussmann et al., 2014, Rall et al., 2010). A decrease in longterm interaction strength may cause an increase of trophic system stability by preventing predator and prey extinctions driven by population fluctuations (Fussmann et al., 2014; Rall et al., 2010; Synodinos et al., 2021). Given that both the functional response parameters and their thermal dependence, and other key parameters such as metabolic rate, have been described to evolve across latitudinal gradients in the longterm (Wang et al., 2021), their rapid evolution during range expansions might be expected (Carbonell et al., 2021; Chuang & Peterson, 2016), and hence also changes in the longterm interaction strength and in the trophic system stability. but not under warming. This provides rare evidence that functional responses can rapidly evolve during range expansions, potentially destabilizing food web dynamics. KEYWORDS ecoevolutionary dynamics, functional response, range expansion, thermal evolution, trophic interactions 13652435, 2022, 12, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2435.14208 by Universidad De Sevilla, Wiley Online Library on [15/06/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 3062 | Functional Ecology CARBONELL et al. The aim of this study is to test for evolutionary changes in the shortterm (functional response) and longterm trophic interaction strengths between an aquatic insect predator (the damselfly Ischnura elegans) and its prey (the water flea Daphnia magna) during the predator's ongoing poleward range expansion in northern Europe. This predator– prey interaction is a textbook example for the functional response in pond food webs (Begon et al., 2006). The prey is native in the predator's expansion region in northern Europe (De Gelas & De Meester, 2005), allowing us to focus on potential evolution of the predator traits (metabolic rate and the functional response parameters) driving changes in the predator– prey interaction strength. We set up a commongarden experiment where predators from core and edge populations at the range expansion front were reared from the egg stage at two temperatures, followed by predation trials to estimate the functional response parameters and by respiration measurements to quantify predator metabolic rate. In a next step, these empirical parameters were used in a predator– prey population dynamic model to assess changes in the longterm interaction strength. The two experimental temperatures (20 and 24°C) reflect the current mean summer water temperature in the study region (20°C) and the future mean temperature assuming 4°C warming by 2100 under IPCC (2021) warming scenario SSP58.5. Based on previous research on the study species (Thompson, 1978; Wang et al., 2021) and another damselfly species (Twardochleb et al., 2020), we expected warming to increase attack rates and shorten handling times thereby increasing shortterm interaction strength, but to decrease longterm interaction strength. However, thermal responses may be more complex and attack rates may increase or decrease depending whether the temperature is below or above the thermal optimum (Englund et al., 2011; Grigaltchik et al., 2012). During this range expansion, larvae from the edge populations evolved a faster development, and with higher food intake and reduced metabolic rate at a higher temperature (Carbonell et al., 2021). Building on this, we here tested for rapid evolution of functional response parameters and their thermal plasticity, and the consequences for the shortand longterm predator– prey interaction strengths during the poleward range expansion. Under the assumption of selection for greater energy needs in the edge populations, we expected the evolution of higher attack rates and/or shorter handling times in the populations from the edge region (based on Urban et al., 2020). 2 | MATERIALS AND METHODS 2.1 | Study system Ischnura elegans (Vander Linden 1820) is a widely distributed damselfly species across Europe with its historical core distribution range going from southern France to southern Sweden (Boudot & Kalkman, 2015). Within its historical core range, I. elegans shows latitudeassociated thermal adaptation in larval life history, behaviour and physiology (e.g. Debecker et al., 2017; Debecker & Stoks, 2019; Dinh Van et al., 2013). During the last few decades, I. elegans is expanding its range northwards in response to global warming (Dudaniec et al., 2018, 2021; Hickling et al., 2005). It has extended its northern range limit by ~300 km in Sweden since 2000 (Lancaster et al., 2015). This would imply the range expansion occurred in ~10 generations, as for Sweden the species seems to have a mixture of a 1year and 2year life cycle at the range front (Szymon Sniegula, pers. comm.). The range expansion in Sweden has been associated with the rapid evolution of an improved cold tolerance, a faster development rate, higher food intake and reduced metabolic rate at high temperatures (Carbonell et al., 2021; Lancaster et al., 2015). This rapid trait evolution is associated with allele frequency changes to facilitate adjustment to the novel thermal conditions at the edge front (Dudaniec et al., 2018). To study the rapid thermal evolution of functional response parameters and the longterm interaction strength in association with the poleward range expansion in Sweden, we selected two regions: a core region in southern Sweden and Denmark, which is within the historical core range (Lancaster et al., 2015, 2016, 2017; Svensson et al., 2005) and the most likely source of the expansion front populations (Dudaniec et al., 2018), and an edge region in central Sweden representing the current poleward expansion front (Figure 1). Both regions are within the historical range of its prey, the water flea D. magna (De Gelas & De Meester, 2005). Based on single nucleotide polymorphisms (including fitnessrelated genes), the populations in the core and edge regions of I. elegans in Sweden are clearly genetically differentiated (figure 2 in Dudaniec et al., 2018). In each region, three replicated populations were sampled in July 2017 (Figure 1, Table S1). We compared the air thermal regimes between the here studied core and edge populations using climatic data from BIOCLIM (Fick & Hijmans, 2017) (see Table S1). In the shallow water bodies, the study species typically occupies water temperatures are highly correlated with air temperatures (e.g. Ali et al., 2016; Prapaiwong & Boyd, 2012). Compared to the studied core populations, mean annual air temperatures (mean ± 1 SE, edge: 5.94 ± 0.29°C vs. core: 7.62 ± 0.35°C), and minimum winter air temperatures (edge: −2.47 ± 0.67°C vs. core: 0.17 ± 0.47°C) are lower in edge populations (as reported in Dudaniec et al., 2018; Lancaster et al., 2015). Mean maximum summer air temperatures (edge: 15.40 ± 0.07°C vs. core: 15.72 ± 0.22°C) and maximum summer air temperatures (edge: 19.20 ± 0.36°C vs. core: 18.90 ± 0.14°C) are similar in both regions. Note that summer temperatures are lower than described in Dudaniec et al. (2018), probably because our study included more coastal populations. 2.2 | Collecting and housing In each of the six populations, 5– 10 mated females of I. elegans (total of 41 females) were collected and placed individually in 100 ml plastic cups with moistened filter paper as oviposition substrate. Filter papers with eggs were sent to the laboratory in Leuven (Belgium) and kept at 22°C for hatching. To enhance initial survival (De Block & Stoks, 2003), during the first 10 days after hatching larvae were 13652435, 2022, 12, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2435.14208 by Universidad De Sevilla, Wiley Online Library on [15/06/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 3063 Functional Ecology CARBONELL et al. placed in groups of 100 in 200 ml plastic containers at the same temperature (22°C) and fed ad libitum with nauplii of Artemia sp. 2.3 | Experimental design To investigate the effects of region (core and edge region) and temperature on life history, physiological and predation traits, a full factorial experiment was set up with three replicated populations per region tested at two rearing temperatures (20 and 24°C). We started the rearing experiment with 100 larvae per combination of region and temperature (total of 400 larvae); the larvae per temperature were as equally as possible distributed among the females per population and the three populations per region. Larvae were reared individually and their source population was tracked during the trials. During the main growth period in spring and summer, larvae in the core region (Southern Sweden and Denmark) typically encounter a mean water temperature of 20°C, while temperatures of 24°C are also encountered, yet in rare events (Debecker & Stoks, 2019). Based on the BIOCLIM data, these thermal patterns are similar in the edge regions. The 4°C temperature difference also matches the predicted global average surface temperature increase by 2100 as predicted by Intergovernmental Panel on Climate Change (IPCC, 2021) scenario SSP58.5. When larvae were 10 days old they were divided in incubators set at 20 and 24°C and a 14:10 L:D photoperiod. Each larva was reared individually in a circular plastic cup (5 cm height, 6 cm diameter) filled with 100 ml dechlorinated tap water. Cups were refilled twice a week to keep the water level constant. Larvae were fed ad libitum with nauplii of Artemia sp. (224 ± 8 nauplii per food ration, mean ± SE, n = 52 food rations) six days a week. When larvae reached the final instar (F0 stage) they were fed daily to meet their higher energy demands. 2.4 | Response variables Key traits related to life history (survival, development rate and growth rate), physiology (metabolic rate) and the functional response parameters (attack rate and handling time) were quantified. For a given larva, all end points were quantified at its rearing temperature. As we here focus on parameters directly related to predator– prey dynamics, methods and results for life history traits are detailed in Appendix S2. Directly after the growth rate estimation, one subset of larvae was used to estimate the functional response, and another subset to measure metabolic rate. Routine metabolic rate (MO2) was determined using closed respirometry (Carbonell et al., 2021). Metabolic rate was estimated on a subset of 15 larvae per combination of temperature and region (total of 60 larvae). For this, oxygen sensitive spots (OxyDot/O2xyDot; Oxysense) were mounted inside 160 ml clear, airtight glass jars. Each jar was supplied with aerated, dechlorinated tap water that was stored in disinfected plastic FIGURE 1 Study populations from the core (red circles) and edge (dark blue triangles) regions of the damselfly Ischnura elegans along the poleward range expansion in Fennoscandia. 13652435, 2022, 12, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2435.14208 by Universidad De Sevilla, Wiley Online Library on [15/06/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 3064 | Functional Ecology CARBONELL et al. buckets in the incubators 12 h prior to the experiment to allow the water to reach the experimental temperature. After sealing the jars, an initial oxygen measurement was taken and glass jars were kept inside an incubator matching the temperature treatment. Oxygen levels in the jars were measured using a Fibox 4 fibre optic oxygen meter (PreSens; Precision Sensing GmbH). After 24 h, a final oxygen measurement was taken to estimate the oxygen consumed by the larvae during that period. To control for background changes in oxygen measurements, three empty glass jars were included in each trial. MO2 was expressed as nmol O2 h−1 STP (standard temperature and pressure) per unit larval wet mass (mg) (Calosi et al., 2013; Carbonell et al., 2017, 2021). To study the predator– prey interaction strength we needed estimates of the attack rates a and the handling times h. Therefore, we reconstructed functional response curves for each of the four combinations of region and rearing temperature. For each functional response curve, we ran predation trials at seven prey densities: 5, 10, 15, 30, 50, 90 and 120 Daphnia/L, based on previous work on the studied species (Thompson, 1978; Wang et al., 2021) and pilot trials that identified prey densities at which food always remained at the end of the trials at the two test temperatures. We used a single D. magna clone from Belgium (Langerodevijver, Huldenberg) in all trials, which was maintained at 22°C. This allowed strong standardization of the prey across all treatments. The choice of a Belgian Daphnia clone can be motivated by the fact that predation rates of Swedish damselfly larvae do not differ when feeding on Daphnia from Belgium or Sweden (De Block et al., 2013). On the other hand, this approach ignored possible differences in thermal adaptation between Daphnia from the damselfly's core and edge regions. Given that we ran predation trials at simulated summer temperatures that were very similar between both regions (Table S1) any such local thermal adaptation of the Daphnia is, however, not expected. Note also that both regions are in the core area of the Daphnia (De Gelas & De Meester, 2005) so no evolution of Daphnia linked to range expansion is to be expected. At each prey density, 8 replicates were run, resulting in a total of 224 predation trials. Each damselfly larva and Daphnia was used in only one predation trial. Prior to predation trials, Daphnia were placed in the arena at the temperature treatment for 2 h (based on Sentis et al., 2016; Thompson, 1978). Trials were conducted in white 2 L (18.0 cm × 13.3 cm × 12.6 cm) plastic containers filled with 1 L dechlorinated water at the treatment temperature. After 24 h of starvation to standardize hunger levels, one damselfly larva was gently transferred to an arena containing one of the seven densities of juvenile (4day old) D. magna. The predation trials directly started after transferring the predator (as in e.g. Sentis et al., 2012, 2016; Wang et al., 2021). Given damselfly larvae immediately start hunting in the arena, any transient stress caused by this transfer is unlikely to affect the functional response parameters. After 24 h, damselfly larvae were removed from the arena and the number of remaining prey counted. For each day of the functional response trials, one control container with 90 daphnids was set at the trial temperature to test natural prey mortality. We chose as control one of the highest densities used in the trials as we did not expect low densities to compromise prey survival. No Daphnia died in these controls. Predation rates were expressed as the number of Daphnia consumed per damselfly larva per 24 h. 2.5 | Statistical analyses and modelling The effects of rearing temperature, region and their interaction on the life history traits and the metabolic rates were tested using linear mixed models (LMMs). Population nested within region was added as a random factor in the models, but it was removed from the final models as it was never significant (pvalue >0.05). Significant Temperature × Region interactions were further explored using Fisher's post hoc LSD t tests with Bonferroni correction; the post hoc results are visualized in the figures. To analyse the functional response curves, the procedure developed by Pritchard et al. (2017) was followed. The fit of type II and type III functional response models was compared based on Akaike's information criterion (AIC). To estimate the attack rate a (initial slope of the functional response; litres/day) and handling time h (inverse of the asymptote of the functional response; days), the type II Rogers's (1972) random predator equation (function frair_fit), which accounts for prey depletion during the predation trials, was fitted as: where N is the number of prey eaten, N0 is the initial prey density, and T is the total duration of the predation trial (in days). The fitted functional response curves were bootstrapped (n = 999) to visualize 95% confidence intervals using the function frair_boot. We pairwise compared the attack rates and handling times between treatment combinations using the ‘indicator variable’ approach (function frair_compare) (e.g. Siepielski et al., 2020). We performed multiple pairwise comparisons between the four treatment combinations, and adjusted pvalues with the FDR method. Attack rates and handling times were used to calculate a new integrated measure of shortterm interaction strength (the per capita impact of the predators on the prey) by combining them into the functional response ratio (FRR) as suggested by Cuthbert et al. (2019): A large attack rate a combined with a short handling time h results in a high FRR value (hence indicating a large per capita effect), while a low attack rate a and a long handling time h gives a low FRR value (a low per capita effect). Following Synodinos et al. (2021), longterm stability (S) was estimated as: (1) N =N0 ( 1−expa[hN−T] ), (2) FRR = a ∕ h . (3) S =− ( 𝜅 − 𝜌 −1) 𝜌−1. 13652435, 2022, 12, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2435.14208 by Universidad De Sevilla, Wiley Online Library on [15/06/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 3065 Functional Ecology CARBONELL et al. S > 0 correspond to a stable equilibrium and S < 0 to oscillations. Kappa, κ, is the ratio of the prey equilibrium density (carrying capacity) without predators to the prey equilibrium density with predators and is thus a measure of longterm interaction strength (see Appendix S3 for more details). Rho, ρ, is the maximal energetic efficiency which corresponds to the ratio between energy gain through food assimilation and energy lost through metabolism assuming a maximum feeding rate (see Appendix S3 for more details). We used our empirical estimates of metabolic rate, handling time and attack rate to calculate k, ρ and S and propagated the errors associated with the estimates of each experimental parameter to obtain their 84% confidence intervals (CIs) (see Appendix S3 for more details). Unlike 95% CIs, a pattern of nonoverlapping 84% CIs matches a difference between two values based on z test at the pvalue <0.05 level (MacGregorFors & Payton, 2013; Payton et al., 2003). Treatment combinations with nonoverlapping 84% CI (analogous to a test for a difference with pvalue <0.05) were considered as significantly different (see e.g. Harris et al., 2020). Statistical analyses and functional response variables and modelling were conducted using R version 3.5.3 for Windows (R Core Team, 2015) using the packages ‘lme4’, ‘lmerTest’, ‘bbmle’, ‘emdbook’, ‘deSolve’, ‘AICcmodavg’, ‘frair’ and ‘propagate’. 3 | RESULTS 3.1 | Predator metabolic rate Predator metabolic rates were not significantly affected by region or temperature (Region: Wald χ2 = 2.0502, df = 1, pvalue = 0.8310; Temperature: Wald χ2 = 0.9884, df = 1, pvalue = 0.6540; Region × Temperature: Wald χ2 = 0.1191, df = 1, pvalue = 0.7310) (Figure 2). 3.2 | Functional response parameters and shortterm interaction strength For both core and edge populations, prey consumption increased with prey density and then reached a plateau following a type II functional response curve at each rearing temperature (Figure 3). The attack rate (a) only significantly decreased at 24°C compared to 20°C in predators from the edge region (Table 1, Figure 4a). Attack rates of predators from core and edge regions were not significantly different at each rearing temperature (Table 1, Figure 4a). The handling time (h) of predators from the edge region was shorter at 24°C compared to 20°C, while that of predators from the core region showed the opposite pattern (Table 1, Figure 4b). The handling time of predators from the core region was shorter than that of predators from the edge region at 20°C but did not differ at 24°C (Table 1, Figure 4a). Based on patterns of overlap of 84% CI intervals, our integrated estimate of shortterm interaction strength, the functional response ratio (FRR), was lower at 24°C than at 20°C for predators from the core region, and was at 20°C lower for the predators from the edge region than those from the core region (Figure 4c). The FRR did not show significant differences between temperatures for predators from the edge region, and between predators from edge and core regions at 24°C (Figure 4c). 3.3 | Longterm interaction strength and predator– prey system stability Based on the overlap of the 84% CI we could not detect significant differences in the longterm interaction strength (κ) nor in the maximal energetic efficiency (ρ) between regions and rearing temperatures. We note, however, that the estimated means of the maximal energetic efficiency were at 20°C about twice as large in core compared to edge populations. All four combinations of region and temperature yielded instable predator– prey longterm dynamics (S < 0) (Figure 5c). At 20°C, the instability was more pronounced for the predator– prey system from the edge region than that from the core region, while at 24°C the instability did not differ between regions. Predator– prey longterm dynamics for the edge region was less stable at 20°C than at 24°C, while differences between temperatures were not found for the core region (Figure 5c). 4 | DISCUSSION To understand and predict the success and ecological impact of species range expansions, we need to consider not only the evolution of species traits associated with range shifts, but also associated changes in trophic interactions with other species and the effects of temperature on these interactions (Diamond, 2018; Jones et al., 2020; Nadeau FIGURE 2 Thermal patterns in metabolic rate (oxygen consumption) of larvae from the core edge regions of the damselfly Ischnura elegans. Given are means with 1 SE. 13652435, 2022, 12, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2435.14208 by Universidad De Sevilla, Wiley Online Library on [15/06/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 3066 | Functional Ecology CARBONELL et al. & Urban, 2019). We therefore studied the rapid evolution of functional response parameters and metabolic rate in the context of an ongoing poleward range expansion. While the functional response parameters and their thermal plasticity showed signals of evolution during the range expansion, this was not the case for metabolic rate. Our results suggest rapid evolution of a lower shortterm impact of the predator on its prey at 20°C to be driven by longer handling times. Instead, the projected longterm interaction strength (estimated as κ) did not seem to change during the range expansion. Nevertheless, the rapid evolutionary changes during the range expansion at 20°C resulted in a less stable predator– prey system from the edge region compared to the predator– prey system from the core region. Because we worked FIGURE 3 Type II functional response curves of larvae of Ischnura elegans from core (a, b) and edge (c, d) regions preying on Daphnia magna at the two rearing temperatures. Delta AICvalues (relative to type II model) refer to the comparison of functional response models II and III. Circles in different colours represent damselfly larvae from different populations of the same region. Shaded areas represent 95% confidence bands. Treatment 1 (T1) Treatment 2 (T2) D Estimate (T1 − T2) SE z pvalue C - 2 0 ° C C - 2 4 ° C a1.58855 0.99270 1.6002 0.1095 h−0.0259 0.00299 −8.6809 <0.0001 E - 2 0 ° C E - 2 4 ° C a3.78866 1.56431 2.4219 0.01544 h0.02833 0.00458 6.1823 <0.0001 C - 2 0 ° C E - 2 0 ° C a−1.27098 1.66572 −0.763 0.4455 h−0.04953 0.00396 −12.521 <0.0001 C - 2 4 ° C E - 2 4 ° C a0.92681 0.80784 1.1473 0.2513 h0.00477 0.00378 1.2605 0.2075 TABLE 1 Pairwise comparisons between functional response parameter estimates for Ischnura elegans damselfly larvae between regions (C = core, E = edge) and rearing temperatures (20 and 24°C) using the ‘indicator variable’ approach (function frair_compare from the frair package). D = difference between the attack rates (a) and handling times (h) between two treatments. pvalues were false discovery rate (FDR)- corrected based on the number of pairwise comparisons (2 comparisons within each region and 2 comparisons between region). Significant values are highlighted in bold 13652435, 2022, 12, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2435.14208 by Universidad De Sevilla, Wiley Online Library on [15/06/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 3067 Functional Ecology CARBONELL et al. with offspring from fieldcollected mothers, we cannot fully exclude a contribution of maternal effects to these changes. Larvae from both regions do, however, show genomic signatures of local adaptation (Dudaniec et al., 2018), indicating rapid evolution to play a role during this range expansion. 4.1 | Patterns in life history traits Consistent with a previous study on this range expansion, we did not observe differentiation in life history between the populations from the core and edge regions at the here studied temperatures of 20 and 24°C (Carbonell et al., 2021). While rapid evolution of life history has been described at moving range fronts (Phillips, 2009), this is not always observed (e.g. Courant et al., 2019). Recent theory predicts that both a faster and a slower life history can evolve at the range front, with the latter being more likely under stressful conditions at the range front (Fitt et al., 2019). Survival and development rates increased when going from 20 to 24°C for both regions (Appendix S1) with 20°C being the mean summer water temperature in ponds in the core region (Debecker & Stoks, 2019). This matches the general pattern that thermal optima are often higher than mean environmental temperatures, probably to avoid being exposed to temperatures above optimal, where performance drops drastically (Deutsch et al., 2008; Martin & Huey, 2008). On the other hand, growth rates did not differ between temperatures, suggesting strong developmental acclimation resulting in socalled thermal compensation (Seebacher et al., 2015), and confirming a previous study on this range expansion (Carbonell et al., 2021). 4.2 | Patterns in the functional response For both regions and at each temperature, prey consumption increased with prey density and reached a plateau following a type II FIGURE 4 Thermal patterns in the functional response parameters of Ischnura elegans from core and edge regions preying on Daphnia magna: (a) attack rate (a) (L day−1), (b) handling time (h) (days), and (c) functional response ratio (FRR = a/h). Given are means with 1 SE for a and h, and with 84% confidence intervals for FRR. Asterisks denote significant differences (FDRcorrected. *<0.05, **<0.01, ***<0.001). For FFR, significant differences (p < 0.05) are based on overlap patterns of the 84% confidence intervals. 13652435, 2022, 12, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2435.14208 by Universidad De Sevilla, Wiley Online Library on [15/06/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 3068 | Functional Ecology CARBONELL et al. functional response curve, thereby confirming previous studies on this predator– prey couple (e.g. Thompson, 1978; Wang et al., 2021). Consistent with the general pattern for ectotherm predators, both functional response parameters were significantly affected by temperature (Haubrock et al., 2020; Sentis et al., 2012; Thompson, 1978; Twardochleb et al., 2020; Wang et al., 2021). More importantly, thermal plasticity for these parameters differed between edge and core populations. Attack rates decreased at the higher temperature in the populations from the edge region but not in those from the core region. Decreased attack rates at higher temperatures have been previously described (Englund et al., 2011; Grigaltchik et al., 2012; Wang et al., 2021). This response suggests that under 4°C warming, larvae from the edge region were less efficient at searching for prey at low densities. Our results relate to those in another damselflyDaphnia predator– prey system where the damselfly Enallagma annexum decreased its capture success above the mean environmental temperature (Twardochleb et al., 2020). According to a metaanalysis by Englund et al. (2011), attack rates first increase with increasing temperatures when the metabolic rate is high enough to maintain energetic demands related with foraging, but then decrease at higher temperatures probably because foraging above the optimal temperature is too costly (Sentis et al., 2012). In addition, temperature effects on attack rates may depend on temperature effects on the Daphnia prey as attack rates depend on predator– prey relative velocities (Novich et al., 2014; Öhlund et al., 2015). However, any such effect on the prey may not explain why the effect of warming on attack rates was only detected in interactions with predators from the edge region. The thermal plasticity of handling times was opposite in populations from the core and edge regions. Populations from the edge region decreased their handling time under warming, indicating 4°C warming to result in higher predation rates at high prey densities. This may suggest that at higher prey densities predators from the edge region switched to a more active foraging style. In contrast, FIGURE 5 Thermal patterns in (a) longterm interaction strength (kappa) and (b) maximal energetic efficiency (rho), and (c) the resulting longterm system stability for the predator– prey systems from edge and core regions composed of larvae of the damselfly Ischnura elegans and the water flea Daphnia magna. Values given are estimated means with 84% confidence intervals. Asterisk denote significant differences (<0.05) based on overlap patterns of the 84% confidence intervals. 13652435, 2022, 12, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2435.14208 by Universidad De Sevilla, Wiley Online Library on [15/06/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License