Introduction

Biological invasions are a major driver of global biodiversity loss through ecological impacts that disrupt the functioning of natural systems (Vilà et al. 2011; Simberloff et al. 2013; Dick et al. 2017). Whilst key drivers of change, the success and impacts of invasive alien species can be mediated by other human-induced alterations, such as habitat and climate changes (Didham et al. 2005; Pyšek et al. 2020). Accordingly, context-dependencies present a major challenge to ecological impact prediction (Ricciardi et al. 2013), with myriad abiotic and biotic factors heightening and dampening invader impacts via emergent effects (Dickey et al. 2020). This is particularly the case in nearshore aquatic ecosystems, where temperature regime shifts can decouple keystone trophic interactions in favour of invasive consumers (Morón Lugo et al. 2020). In highly productive coastal waters in particular, changes to temperature, salinity, pH, dissolved oxygen and eutrophication can occur rapidly, potentially disrupting ecosystem structure and function (Boyd et al. 2018). These systems are also especially prone to invasion, with increasing intensity of globalised trade and transport networks connecting discrete coastal habitats and their species communities (Casties et al. 2016).

Invasive predators can be particularly damaging to the stability of native assemblages through an array of trophic (i.e., density-mediated) and non-trophic (i.e., trait-mediated) interactions (Platvoet et al. 2009; Anton et al. 2020). Advances have recently paralleled quantifications of invasive species functional responses with the magnitude of their in-field ecological impact (Dick et al. 2014). The functional response approach quantifies resource use of a consumer as a function of resource density, allowing for comparative elucidation of density-dependences in ecological impacts under a range of abiotic variables (e.g., warming; Wasserman et al. 2018) and species compositions (e.g., presence of multiple predators; Sentis and Boukal 2018). Classically, the functional response has been characterised discretely into three forms: types I (linear), II (hyperbolic) and III (sigmoidal), whereby type II feeding relationships can impart greatest impact due to high feeding rates at low resource (e.g. prey) densities (Holling 1959; Hassell 1978). However, impact predictions for invasive alien species have often been grounded in comparisons of single predators, which ignore emergent multiple predator effects (Sih et al. 1998; Médoc et al. 2013), with predators rarely present singularly in natural communities. For invaders, predator–predator interactions might thus profoundly mediate their impacts through emergent antagonisms and synergisms that can dampen or exacerbate ecological impacts, respectively (Wasserman et al. 2016). In turn, the strengths of multiple predator interactions might be further altered by temperature regime and prey density (Sentis et al. 2017), given the primacies of temperature and prey availability in foraging behaviour (Englund et al. 2011), necessitating quantitative investigations into emergent effects among these parameters.

The present study thus quantifies the influence of warming and predator density on the functional response of a widespread invasive predator that has invaded marine, brackish and freshwater environments in recent decades i.e. Gammarus tigrinus (Crustacea: Amphipoda) Sexton, 1939. Native to North America, this thermophilous and euryhaline species has invaded many European aquatic ecosystems (Casties et al. 2016; Paiva et al. 2018). Indeed, this species is one of few invasive gammarids that has successfully transitioned between marine and freshwater environments (Cuthbert et al. 2020), and is widely distributed and spreading in the brackish waters of the Baltic Sea (Rewicz et al. 2019). Gammarus tigrinus is known to actively feed on a range of aquatic invertebrates and prioritises carnivory to maximise energy intake, irrespective of temperature (Pellan et al. 2016). For this experiment, we used larval chironomid prey. These sediment-dwelling, free-swimming prey occupy benthic areas of waterbodies, and are a widespread and diverse group in the Baltic Sea—estimated to comprise approximately 30% of the macrozoobenthos species in the system, with as many as 230 species there (Ojaveer et al. 2010). Given that temperature regime is closely aligned with feeding rates (Englund et al. 2011), and gammarids within the same guild are known to interact (Platvoet et al. 2009), we expected (1) warming to increase gammarid functional responses towards live prey, and (2) intraspecific multiple predator effects to mediate ecological impacts of invasive gammarids differentially between temperatures, given increasing foraging intensities with warming that exacerbate predator–predator interactions.

Materials and methods

Animal collection and experimental design

The predators, G. tigrinus, were obtained from Travemünde, Lübeck, Germany (53°83′ N 10°64′ E) during summer 2017 and transported to a controlled environment chamber (18 °C; 12:12 light and dark regime) at GEOMAR Helmholtz Centre for Ocean Research Kiel, Germany in source water. The salinity at the time of sampling Travemünde was 10 ppt, owing to adjacent freshwater inputs from an inland waterway that reduce the salinity of the site. A kick-net was used to collect gammarids from within macroalgae at nearshore parts of the site. In the chamber, gammarids were housed in 56 L aquaria (~ 100 individuals per aquaria), after identification, containing filtered Baltic Sea water from Kiel Fjord that was mixed with tap water to reach the desired ambient salinity of the collected population (10 ppt). Water from the Kiel Fjord was used for practical purposes over the acclimation period, however, laboratory conditions were calibrated to reflect those of the sampling location as much as possible. In the aquaria, water was continually pumped through a filtration system in each tank, and each tank contained a mixture of artificial floating vegetation and stones to provide habitat structure. Gammarids were housed in laboratory conditions and acclimated for multiple generations (> 1 year) prior to use. The prey, live larvae of chironomids, were obtained commercially (ZOO & Co. Knutzen, Kiel). Gammarids were fed a mixture of ground commercial fish and shrimp flakes, supplemented with chironomid larvae in the weeks prior to experimentation. Gammarus tigrinus has previously been shown to readily consume chironomid larvae (Pellan et al. 2016).

In the controlled environment chamber, functional responses of G. tigrinus (total length: 0.9–1.1 cm) towards live chironomid larval prey (total length: 0.6–1.1 cm) were quantified factorially under two temperatures (18 °C, 24 °C), five prey densities (2, 4, 8, 16, 32) and in the presence of single or double predator treatments (1 or 2 size-matched gammarid individuals). Controls consisted of prey in the absence of predators under each temperature and density treatment, to quantify non-predatory background mortality rates. At least three replicates were conducted per experimental group; all treatments were initially replicated four times, however, in some cases predator mortality or moulting during the experiment required us to exclude these replicates. Gammarids were unfed for 24 h prior to experimentation to standardise hunger levels in 5 L aquaria containing artificial stones for refuge, separate to the experimental aquaria.

Experiments were conducted in 500 mL plastic arenas of 12 cm diameter containing 10 ppt water, within water baths maintained at the nominal experimental temperatures (18 °C or 24 °C). This temperature range is within the magnitude of daily fluctuations of temperature recorded in the study area (Pansch and Hiebenthal 2019), with the lower temperature reflecting ambient conditions in the Kiel Fjord (Morón Lugo et al. 2020). The higher temperature reflects changes within nearshore areas of the Baltic Sea that occur sporadically throughout the year, based on empirical data from the last two decades, owing to heatwaves which are becoming increasingly frequent (Pansch et al. 2018). Both predators and prey were acclimated separately in the 500 mL arenas to the higher temperature for 2 h at a rate of 1.5 °C 30 min−1. Following acclimation, predators were added to each prey density and allowed to feed for 24 h (i.e. the duration of the experiment). Owing to laboratory logistics, feeding trials were run on 7 separate days, with replicates randomised within each temperature grouping per day. After this period, gammarids were removed and remaining live prey enumerated to quantify those killed via predation. Dissolved oxygen levels were monitored in both temperature treatments, and these did not fall below 80% saturation in either case. We did not observe any evidence for predation between gammarids in intraspecific multiple predator treatments.

Statistics

Median levels of prey mortality in predator-free controls at each prey density and temperature were used to account for natural prey mortality in the calculation of consumption rates in those treatments containing predators (i.e., via subtraction of control mortality from predator-driven mortality). Counts of prey killed were analysed using generalised linear models assuming a Poisson error distribution and log link, as a function of temperature, prey density and predator numbers. An information theoretic approach was used to select the model which minimised information loss via second-order derivations of Akaike’s information criterion (AICc; Burnham and Anderson 2002). A selection table was compiled with all possible combinations of explanatory terms in the model (Barton 2020), with models then ranked based on AICc, whereby lower AICc indicates a better fit. Models with ΔAICc ≤ 2 were considered interchangeable. The final, selected model thus comprised that with the lowest AICc, with analysis of deviance with type III sums of squares then used to compute coefficients of the best fitting model (Fox and Weisberg 2019). Tukey comparisons via estimated marginal means were used for pairwise testing where a term was significant (Lenth 2020).

Gammarid functional response types (I, II or III) were categorised at the single-predator density using binomial generalised linear models with logit links for each temperature treatment separately, with consumption rates (proportion of prey killed) analysed as a function of initial prey density. Because prey were not replaced as they were killed, we fit Rogers’ random predator equation (Rogers 1972; Pritchard et al. 2017) to the data to model functional response parameters (attack rate and handling time), assuming the data were properly described by a type II functional response curve (Haddaway et al. 2012; Rosenbaum and Rall 2018):

$$N_{{\text{e}}} = N_{0} \left( {1 - {\text{exp}}\left( {a\left( {N_{{\text{e}}} h - T} \right)} \right)} \right),$$
(1)

where Ne is the number of prey eaten, N0 is the initial density of prey, a is the attack rate, h is the handling time and T is the total experimental period. Because Ne appears on both sides of Eq. 1, the solution was found using Lambert’s transcendental equation (Bolker 2008). The difference (△) method was then used to compare functional response parameters (a and h) pairwise between temperatures (Pritchard et al. 2017).

To quantify the strength of non-trophic interactions in treatments with multiple gammarid predators between both temperatures, we first quantified the strength of trophic interactions that included both tropic (i.e., predatory) and non-trophic (i.e., trait-mediated) effects from experimental observations:

$${\text{IS}}\left( {P,Z} \right) = \frac{{N_{{\text{P}}} - N_{{\text{P,Z}}} }}{{N_{{\text{P}}} }},$$
(2)

where NP and NP,Z are the proportions of live prey at the beginning and end of the experiment, respectively. Second, using single-predator functional response parameters (a and h) quantified as per Eq. 1, we fit a population-dynamic model to simulate expected predation rates in the absence of non-trophic interactions (i.e., without intraspecific multiple predator effects) (McCoy et al. 2012; Sentis and Boukal 2018):

$$\frac{{{\text{d}}N}}{{{\text{d}}t}} = - \mathop \sum \limits_{i = 1}^{n} f_{i} \left( N \right)P_{i} ,$$
(3)

where N is the prey population density, Pi (i = 1, 2,…, n) are the population densities of predators i and fi(N) is the functional response of predator i (i.e., Eq. 1). The equation was integrated over time to obtain expected numbers of surviving prey at each temperature treatment and prey density using a Latin hypercube sampling algorithm (Soetaert et al. 2010; Soetaert and Petzoldt 2010). Last, to quantify the strength of non-trophic interactions, the simulations from Eq. 3 (i.e., excluding non-trophic interactions) were subtracted from the observations from Eq. 2 (i.e., including non-trophic interactions) to quantify intraspecific multiple predator effects. Here, negative resulting values would indicate antagonistic predator–predator interactions given that predictions exceed observations. Conversely, synergisms would be indicated by positive values (i.e., observations exceeded predictions). Non-trophic interactions were analysed using a Scheirer–Ray–Hare test as a function of temperature treatment and prey density. All analyses were computed in R v4.0.2 (R Core Team 2020), with significance inferred at an α of 0.05.

Results

The best fitting model of prey consumption included predator density, temperature and prey density as single and interacting terms, but excluded the three-way interaction among these variables. Intraspecific multiple predators consumed significantly more than predators present singularly (Table 1). Predator density and temperature effects, however, significantly interacted to influence prey consumption (Table 1); warming had a significant, positive effect considering single (p = 0.001) but not multiple (p = 0.821) predators. The effects of temperature were therefore different according to predator density, becoming less apparent as predator density increased. Prey density was significant in both interaction terms with predator density and temperature (Table 1). All other single terms were not statistically clear (Table 1).

Table 1 Analysis of deviance table with type III sums of squares considering Poisson generalised linear model of numbers of prey consumed as a function of predator density (Predator), temperature (Temperature) and prey density (Density)
Full size table

Across both temperatures, G. tigrinus consumption rates related negatively with increasing prey density (Table 2), and thus hyperbolic type II functional responses were exhibited (Fig. 1). Attack rates tended to decrease with temperature, whilst handling times tended to shorten and thus maximum feeding rates increase (Table 2). Under single-predator densities, attack rates did not significantly differ between temperatures, owing to high standard error, whilst handling times were significantly shorter with warming (Table 3; Fig. 1). Accordingly, functional response initial slopes were similar between temperatures at low prey densities, yet more marked differences in magnitudes were found at high prey densities (Fig. 1).

Table 2 First-order terms from binomial generalised linear models considering prey killed as a function of prey density across temperature and predator treatments by Gammarus tigrinus. Functional response attack rates, handling times and maximum feeding rates from the random predator equation
Full size table
Fig. 1
figure 1

Functional responses of single Gammarus tigrinus under 18 °C and 24 °C temperatures. Shaded areas are non-parametric bootstrapped (n = 2000) 95% confidence intervals. Points are raw data

Full size image
Table 3 Results from difference (Δ) method comparing functional response parameters of single predators between temperatures, alongside standard errors (SE)
Full size table

The strength of non-trophic interactions between intraspecific multiple predators was affected by a significant interaction between temperature and prey density (H = 17.278, p = 0.002). At low prey densities (2 and 4 ind. arena−1), consumption rates were well-predicted from simulations, and thus no evidence for non-trophic interactions was discerned (Fig. 2). At intermediate densities (8 ind. arena−1), non-trophic interaction strengths were negative under both temperature treatments, indicating interference effects among predators. At highest prey densities (16 and 32 ind. arena−1), non-trophic interactions were, in turn, dependent on temperature, with positive values (i.e., synergisms) observed at 18 °C and negative values at 24 °C (i.e., antagonisms). Thus, the nature of intraspecific multiple predator effects was mediated by emergent interactions between temperature and prey density, with warming driving interference effects at high prey densities.

Fig. 2
figure 2

Non-trophic interaction strengths of multiple Gammarus tigrinus across prey densities and between warming treatments. The solid horizontal line indicates expected values in the absence of non-trophic effects. Positive values indicate synergisms, and negative values antagonisms, between multiple predators

Full size image

Discussion

Ecological impacts from a notorious invasive alien species in marine, brackish and freshwater habitats were shown to be mediated through intraspecific multiple predator effects between two temperatures and across a prey density gradient in the present study. Whilst warming generally intensified interaction strengths with a representative benthic prey (Pellan et al. 2016), corroborating previous studies on gammarid congenerics (Laverty et al. 2017), emergent effects that mediate invader impacts were found according to temperature and prey density. In particular, intraspecific multiple predator interactions eroded the magnitude of positive temperature effects on feeding rates. As such, despite increasing numbers of invasive gammarids heightening overall predation rates, simulations based on individual predator feeding indicated predator–predator antagonisms at higher temperatures and synergisms at low temperatures as prey densities increased. Empirically, given that functional responses are predictive of known ecological impacts (Dick et al. 2017), these results indicate that warming will worsen invader effects, but density-dependent predator–predator interacions could mitigate these impacts. Nevertheless, consistently efficient resource use by G. tigrinus at low prey densities may result in a stronger likelihood of extirpation for rare prey populations (Dick et al. 2014). Given that temperature oscillations can be of a magnitude of 21 °C in the study area (Western Baltic Sea) between winter and summer (Pansch et al. 2018), with up to 7 °C changes reported within a single day due to day-night transitions and upwelling events (Pansch and Hiebenthal 2019), our results show invader impacts could change rapidly over short timescales in these and other coastal regions.

The strength of non-trophic interactions (i.e. multiple predator effects) was shown to differ across prey densities in the present study under both temperature regimes tested. These findings corroborate previous studies that have identified such prey density-dependencies (Sentis et al. 2017), and may be an artefact of changeable predator–predator interactions as resource availability shifts. Specifically, at low prey densities, all prey are typically extirpated in non-prey replacement experimental designs, and thus there is little capacity to detect non-trophic interactions. At intermediate prey densities, competition between predators for limited resources is high, resulting in antagonisms that were detected at both temperatures in the present study. Conversely, at high prey densities, prey are abundant and thus not extirpated, with predator–predator interactions and thus multiple predator effects potentially less pertinent (Sentis et al. 2017). However, our study found temperature to additionally mediate these effects at high prey densities, whereby predators interacted positively at lower temperatures and negatively at higher temperatures. We posit that this is intertwined with predator activity rates, with foraging activity generally increasing in gammarids at higher temperatures (Laverty et al. 2017), potentially increasing predator–predator encounters and interference effects. Indeed, other notorious marine bioinvaders have similarly been shown to exhibit density-dependences in intraspecific multiple predator effects, such as lionfish (De Roy et al. 2020), yet these trends might vary according to foraging traits of individual invasive species (Wasserman et al. 2016). The effects of temperature on encounter rates thus require further investigation.

Invasions in aquatic systems by invasive alien gammarids have exhibited unidirectional patterning, with the majority of species moving from the brackish waters of the Ponto-Caspian region to Eurasian and North American freshwaters globally (Cuthbert et al. 2020). Gammarus tigrinus, however, is an exception to this phenomenon, having successfully invaded marine, brackish and freshwater ecosystems from its native range, that also spans salinities from freshwater to marine (Paiva et al. 2018). The wide-ranging salinity and temperature tolerance of this species thus heightens its invasion success, with this gammarid continuing to spread through the Baltic Sea (Rewicz et al. 2019). In turn, it is likely to impart greater impacts on prey communities in future as coastal waters warm, heatwaves intensify, and coastal regions become more interconnected globally. Whilst our results provide novel insights into the implications of warming and biotic contexts for invader ecological impacts, further work is required that incorporates additional predator densities, taxonomic groups and environmental change scenarios, as well as invader reproductive and prey responses to climate change.