Forecasting impact of existing and emerging invasive gobiids under temperature change using comparative functional responses

Round goby Neogobius melanostomus (Pallas, 1814) and western tubenose goby Proterorhinus semilunaris (Heckel, 1837) have recently expanded substantially beyond their native ranges, posing a threat to freshwater and brackish ecosystems. Both species exert a detrimental effect on fish community structure directly via predation on fish eggs and indirectly via alterations of food webs. While the impact of N. melanostomus is widely documented, P. semilunaris receives little attention and its effect on aquatic communities has not yet been quantified. We compared aspects of their predation on common carp Cyprinus carpio Linnaeus, 1758 larvae at 20 °C and 25 °C using the functional response (FR) approach, which has been developed and applied to forecast invader impact. Intra-specific comparison showed no significant temperature dependency on assessed FR parameters, attack rate and handling time, suggesting broad temperature tolerance of both tested predators. Proterorhinus semilunaris demonstrated a marginally higher attack rate at 20 °C compared to N. melanostomus. The handling times did not differ between predators. Proterorhinus semilunaris showed a lower maximum feeding rate at 25 °C compared to the rest of treatments suggesting lower temperature optima compared to N. melanostomus. Both predators showed substantial maximum feeding rates, which may impact recruitment of native fish species.


Introduction
With continuing growth of international commerce, travel, and transport, ecosystems are under increasing threat from invasions of non-indigenous organisms, which show no signs of saturation (Hulme 2009; Seebens et al. 2017).Biological invasions seriously impact global biodiversity (Blackburn et al. 2004) and ecosystem processes (Ehrenfeld 2010), causing major social and economic problems (Pimentel et al. 2001;Mazza et al. 2014).Global homogenization of previously separated biota has especially pronounced consequences for aquatic environments (Moorhouse and Macdonald 2015).
Neogobius melanostomus causes major alterations in newly invaded areas.In addition to critical food web disruptions and shifts in trophic levels (Krakowiak and Pennuto 2008;Almqvist et al. 2010), there may be habitat and diet overlap with native ichthyofauna (Skóra and Rzeznik 2001), causing a significant decrease in, or extermination of, native fishes (Janssen and Jude 2001).Impact of P. semilunaris, on the other hand, has not yet been quantified and only potential threats and indirect effects are mentioned (Kocovsky et al. 2011).Both N. melanostomus and P. semilunaris are opportunistic feeders capable of adapting to locally abundant food sources (Brandner et al. 2013;Všetičková et al. 2014).Although low consumption of fish eggs and juveniles by P. semilunaris and N. melanostomus have been reported (French and Jude 2001;Vašek et al. 2014), direct predation by non-indigenous gobies on fish eggs and/or larvae in both laboratory and field conditions have also been observed (Thompson and Simon 2014;Všetičková et al. 2014) indicating these resources can form a substantial part of the diet during specific events (e.g.spawning season).Therefore, it is crucial to understand the interactions of these species with potential prey.
Predation is a key process in the population dynamics of lower trophic levels that influences the entire food web regardless of ecosystem type (Holling 1959;Bax 1998).Resource availability (prey density) is a crucial determinant of feeding rate and has been characterized by the use of the functional response curve (FRs; Solomon 1949;Holling 1959).The form and magnitude of FRs are important aspects of consumer-resource interactions and community dynamics, respectively (Murdoch and Oaten 1975;Juliano 2001).Invasive species often display elevated FRs compared to native or low-impact non-native ecologically analogous species and therefore comparative FR is a useful tool for comparison and assessment of the impact of newly introduced species (Dick et al. 2013;Alexander et al. 2014;Dick et al. 2014;Xu et al. 2016;Laverty et al. 2017).
Temperature is considered an important driver of nearly all biological processes, encapsulating individual and population growth rates as well as development rates and lifespans (Brown et al. 2004).Based on the current progression of climate change, we can expect a substantial alteration in consumerresource relationships worldwide (Clavero and Garcia-Berthou 2005;Pereira et al. 2010;Schulte et al. 2011).Petchey et al. (1999) documented a decline in numbers of top predators and herbivores in warmed ecosystems.Changes in this structure may drive cascading effects at population and community levels (Brown et al. 2004;Petchey et al. 2010).Warming may also substantially alter primary and secondary production in aquatic ecosystems (Richardson and Schoeman 2004).While N. melanostomus has a wide temperature tolerance, ranging from −1 to 30 °C (Moskal'kova 1996), maximum consumption is estimated to be within the range of 23 and 26 °C, after which consumption rates decline abruptly (Lee and Johnson 2005).Studies on P. semilunaris thermal requirements are scarce but O 'Neil (2013) suggests the thermal optima is lower than that of N. melanostomus.
In the present study, FRs were applied to compare and forecast the impacts and foraging capacity of both established and emerging invaders, N. melanostomus and P. semilunaris, upon soft-bodied pelagic prey.The impact of both predators was tested at two temperature levels to explore the relative effects at 20 °C, representing the likely spawning temperature of cyprinid fish prey, and 25 °C, representing the temperature of estimated maximum consumption in N. melanostomus.It was hypothesized that N. melanostomus would exhibit elevated FR at 25 °C compared to P. semilunaris due to higher temperature optima of the former species (Lee and Johnson 2005;O'Neil 2013).

Predator acquisition and acclimatization
Western tubenose goby Proterorhinus semilunaris (Heckel, 1837) and round goby Neogobius melanostomus (Pallas, 1814), were collected using a backpack pulsed-DC electrofishing unit (FEG 1500, EFKO, Leutkirch, Germany) in late April 2016 from recent colonized localities: P. semilunaris from the Jevišovka River (48º49′27″N; 16º27′59″E) and N. melanostomus from the Elbe River (50º39′7″N; 14º2′41″E).Fish were transported to the experimental facility of the Institute of Aquaculture and Protection of Waters, Faculty of Fisheries and Protection of Waters, University of South Bohemia in České Budějovice, Czech Republic.Each of these predator species was distributed between two identical recirculating aquaculture systems (RAS) filled with dechlorinated tap water (20.3 ± 0.1 °C, pH 7.4-7.6,oxygen level > 95% saturation).Size classes of both species were not matched (Table 1).During the 30-day acclimatization period, the temperature was maintained at two levels for both species, 20 °C (20.4 ± 0.3 °C) and 25 °C (24.9 ± 0.3 °C).The latter temperature was reached raising the temperature by 0.1 °C h -1 over the first 2 days of the acclimatization period to minimize stress caused by rapid heat change.The oxygen level was > 95% saturation and pH was 7.0-8.0 in both temperature treatments.Temperature, dissolved oxygen, and pH were measured twice daily during the acclimatization period with an HQ40d digital multi-meter (Hach Lange Eight prey densities (8, 20, 45, 100, 180, 290, 420, 550 larvae per box) were used.Carp larvae were placed into each box 1 h prior to introduction of predators to allow recovery from handling stress.Each combination of temperature, prey density, and predator treatment was replicated six times (192 boxes with predators).Five control replicates, with no predator, at each combination of prey density and temperature were included to assess background mortality of prey (80 control boxes).Predators were individually and randomly released into boxes with prey.After 24 h predators were removed and placed in a marked box to allow matching to preyed larvae, and the remaining prey were counted to derive those killed during the experiment.Each predator was used only once to avoid bias caused by experience.Total and standard lengths of predators, measured by digital caliper to the nearest 0.01 mm, and weight to the nearest 0.1 g (Kern 572-35, Kern and Sohn, Germany) were recorded after the feeding experiment (Table 1).

Statistical analysis
We employed logistic regression to examine predation rate as a function of prey density and decipher the form of the functional response (FR) in each temperature regime: (1) where N e is the number of prey consumed per gram of predator, N 0 is the initial prey density and P 0 , P 1 , P 2 , and P 3 are the intercept, linear, quadratic, and cubic coefficients, respectively, estimated by maximum likelihood (Juliano 2001).If P 1 < 0, the proportion of prey killed declined monotonically with the initial prey density, indicating a type II FR.If P 1 > 0 and P 2 < 0, the proportion of prey killed is a unimodal function of prey density, corresponding to a type III FR (Juliano 2001).We used type II Rogers' random predator equation (Rogers 1972), which accounts for prey depletion, to model all functional responses: where N e is the number of prey consumed per gram of predator, N 0 is the initial prey density, a is predator attack rate classically defined as searching efficiency, h is predator handling time defined as time spent pursuing, subduing, and consuming each prey item plus the time spent preparing to search for the next prey item, and t is the duration of the experiment in days.Owing to the implicit nature of the random predator equation, we used the Lambert W function to solve Equation (2) (for further details see Bolker 2008): We also tested whether predator a and h values were influenced by temperature and whether a and h values differ among predators within a temperature level.For each test we evaluated four FR models covering all possible combinations of temperature dependence and predator comparison in both FR parameters (Sentis et al. 2015) using the maximum likelihood method implemented in the "bbmle" package (Bolker 2008) and selected the most parsimonious model based on the lowest AICc values (Burnham and Anderson 2002).
Weight-specific maximum feeding rate was derived from the FR models (Dick et al. 2013;Alexander et al. 2014).Bootstrapping was used to generate multiple estimates (n = 30) of the parameter of the maximum consumption rate: (4) followed by full factorial comparison (i.e."temperature" and "predator species" as explanatory variables) using GLM (due to non-normal data distribution) with quasi-Poisson error distribution and simplification via step-deletion process.Tukey HSD post-hoc tests were employed to assess significant differences among treatment means.All analyses were implemented in R version 3.2.5 (R Core Team 2016).

Results
In all controls, prey survival was > 98% after 24 h, thus mortality was attributed to predation, and data were not corrected for natural mortality.

Functional response
Based on the results of logistic regression, both predators showed type II FR at both temperatures (Linear parameters were negative, and the proportion of fish larvae eaten decreased with increasing prey density, Figure 1, Supplementary material Table S1).Weight-specific predator intake rate increased with prey density, reaching an asymptote corresponding with the maximum number of prey eaten in 24 h.

Attack rate and handling time
Estimates of FR parameters relative to temperature regime and predator species are provided in Table S2.
Based on the most parsimonious model, we found neither attack rate (a) nor handling time (h) to significantly differ with temperature within either species (Table 2).In inter-specific comparison, a was found to be marginally higher in P. semilunaris at 20 °C (Table 3).Handling time did not differ between species in either temperature regime (Table 3).

Maximum feeding rate
Model selection revealed significant positive interactions among temperature and predators species (Table 4).Weight-specific maximum consumption of N. melanostomus at 20 °C and 25 °C and P. semilunaris at 20 °C showed no significant difference (mean ± SE: 84.5 ± 1.2; 85.1 ± 1.2 and 82.0 ± 1.1 fish larvae per gram of predator, respectively).At 25 °C, however, P. semilunaris maximum feeding rate was significantly lower compared to the other treatments (70.7 ± 1.0 fish larvae per gram of predator; Tukey HSD; Table 5).

Discussion
Many studies have been devoted to the negative impact of N. melanostomus (Kornis et al. 2012), whereas P. semilunaris has received little attention and only indirect effects are mentioned (Kocovsky et al. 2011).Calculation of Functional Response allows effective assessment and comparison of impacts of invader/native and invader/invader pairings (Dick et al. 2014;Xu et al. 2016).The FR defines the relationship between resource availability (prey density) and consumption pattern (Solomon 1949;Holling 1966) and is an essential component of ecosystem modeling (Petchey et al. 2008).In our study, both predators showed type II FR.The shape of the FR curve corroborated previous findings for N. melanostomus (Dubs and Corkum 1996;Fitzsimons et al. 2006) and for other piscine predators (e.g.Alexander et al. 2014); hence type II FR is a common foraging  pattern for this ecological group.Consequently, both predators may have a destabilizing effect on population dynamics of cyprinids through high prey exploitation at low prey densities, typical for type II FR (Murdoch and Oaten 1975).Although the P. semilunaris FR curve showed a steeper gradient, N. melanostomus curves reached slightly higher asymptotes at both temperatures.
Temperature is believed to be one of the most important environmental factors influencing predatorprey interactions (Petchey et al. 1999;Brown et al. 2004).Although published data suggest the maximum consumption rate of N. melanostomus falls within 23-26 °C, alteration from 20 to 26°C produced only a slight increase (Lee and Johnson 2005).This is consistent with our results, where N. melanostomus showed no significant differences between 20 °C and 25 °C in either attack rate or handling time.A similar conclusion is drawn for P. semilunaris.Attack and maximum feeding rates of ectotherms are typically characterized by humpshaped temperature dependency and, while sharing the same temperature optima, attack rate is dynamically affected, both above and below the optimal temperature (Englund et al. 2011).Our results therefore imply that tested temperatures fall within the optima of both predator species, as we observed no significant intraspecific difference in attack rate between 20 °C and 25 °C.On the other hand, South et al. (2017) recorded higher FR curves of invasive lionfish Pterois volitans (Linnaeus, 1758) at 26 °C compared to 22 °C, whilst attack rates were similar.
Invasive species are often characterized by a higher maximum feeding rate, higher attack rate, and shorter handling time compared to native or low-impact non-indigenous species (Xu et al. 2016).These trends are apparent in fishes (e.g.Dubs and Corkum 1996) as well as in other aquatic organisms including amphipods (Bollache et al. 2008), crayfish (Renai and Gherardi 2004), and mysids (Dick et al. 2013).Within invader/native comparisons, Alexander et al. (2014) found similar attack rates in benthic dwelling fish species attributed to similar predatory strategies, whereas in pelagic fishes, the native showed an even higher attack rate.In our study, P. semilunaris showed a marginally higher attack rate at 20 °C compared to N. melanostomus (owing to the slightly higher AICc value), while no difference was observed at 25 °C.In a comparison between non-indigenous N. melanostomus and the analogous native predator European bullhead (Cottus gobio Linnaeus, 1758), Laverty et al. (2017) found higher FR asymptotes in N. melanostomus, indicating lower handling times.In comparing two alien predators, we observed no significant differences in handling time at either tested temperature, suggesting that the per capita impact on native communities posed by both predator species is similar, although little attention has yet been given to the impact of P. semilunaris.However, introducing Relative Impact Potential (RIP), Dick et al. (2017) suggested that both per capita effect (FR) and numerical response (abundance) should be included when estimating invader ecological impact.
In general, P. semilunaris do not reach as high an abundance as N. melanostomus in invaded ecosystems (Dopazo et al. 2008;van Kessel et al. 2016;Janáč et al. 2018), and consequently the former species is no longer considered an important threat to native communities (Borcherding et al. 2011).On the other hand, P. semilunaris reaches high densities in shallow macrophyte-rich habitats (Kocovsky et al. 2011), where it can pose serious threat for phytophilic fauna, including cyprind larvae.Interestingly, N. melanostomus abundance in nearshore areas is positively correlated with temperature up to 24 °C (Dopazo et al. 2008) and hence, the RIP can be higher under such conditions.Compared to an analogous native predator C. gobio, N. melanostomus consumes significantly more prey (Laverty et al. 2017).In laboratory conditions, Ray and Corkum (1997) documented N. melanostomus average consumption of 1.02 g of zebra mussel, Dreissena polymorpha (Pallas, 1771), wet weight per 24 h in larger specimens (SL 70-84 mm), whereas smaller fish (SL 55-69 mm) consumed ~20% that amount.This finding suggests that juveniles may depend on other food sources, possibly fish larvae or other soft-bodied organisms, rather than hard-shelled mollusks.In our study, the weight-specific maximum feeding rate of fish larvae by N. melanostomus  showed no difference between the tested temperature regimes.Similar maximum feeding rate was reached by P. semilunaris (SL 35-58 mm) at 20 °C.However, at 25 °C, P. semilunaris maximum feeding rate was significantly lower compared to the other treatments.These results correspond with the modeled weight-maximum consumption curve of N. melanostomus, in which smaller specimens showed a higher weight-specific consumption rate (Lee and Johnson 2005).The model for maximum consumption rate in P. semilunaris has not been developed, but, based on our results, we can expect a lower temperature optima of this species, corroborating O 'Neil (2013).Based on the seasonal temperature course in a Central European lowland river (River Elbe Board, State Enterprise, unpublished data), the investigated species can potentially reach maximum feeding rate from late May/early June through early October.

Conclusions
Both N. melanostomus and P. semilunaris, readily consumed cyprinid fish larvae in laboratory conditions.No significant differences in handling time relative to temperature were recorded suggesting a similar per capita impact on native communities.Considering the large numbers of fish larvae consumed by both studied species, the pelagic lifestyle of juveniles (Jůza et al. 2015), and abundance in invaded water bodies (van Kessel et al. 2016), the impact on native fish populations by their direct predation may be considerable.The higher maximum feeding rate of N. melanostomus at 25 °C compared to P. semilunaris indicates higher temperature optima of the former species.Further studies are needed to investigate multiple predator effects of N. melanostomus and P. semilunaris, as they are often recorded sympatrically (van Kessel et al. 2016) and thus interact with one another (Johnson et al. 2009;Wasserman et al. 2016).Despite significant differences in predator sizes between treatments, we believe that such minor differences are not biologically relevant.Yet, size-matching of tested species should be considered in future experiments to avoid these biases.

Figure 1 .
Figure 1.Relationship between prey density and number of prey eaten per gram of predator for Neogobius melanostomus and Proterorhinus semilunaris.Data are shown as mean ± SE overlaid by the model prediction (blue solid line = 20 °C, red dashed line = 25 °C).

Table 1 .
Biometric data of Neogobius melanostomus and Proterorhinus semilunaris used in experiments at two temperatures.Total length (TL), standard length (SL), and weight (W).Data are presented as mean ± SE.Different letters in particular columns denote significant differences among groups (Kruskal-Wallis, Multiple comparison of mean ranks) at α = 0.05.
Lee and Johnson 2005) and Peňáz 1994)2000)1758 at the early protopterygiolarval stage (3-4 days post hatching:Balon 1975)were used as prey.This stage was chosen as it represents larva size and behavior of cyprinids (seeOsse and van den Boogaart 1995;Çalta 2000), representatives of the major fish communities in many European reservoirs, as well as regulated lowland rivers(Kubečka 1993;Jurajda and Peňáz 1994).Common carp larvae were obtained from artificial reproduction and incubation at the Genetic Fisheries Center, Faculty of Fisheries and Protection of Waters, University of South Bohemia.weretested,20°C and 25 °C.The former represented the temperature of natural reproduction of common carp(Horváth et al. 1984), while the latter is within the reported range in which N. melanostomus exhibits maximum consumption rate (23-26 °C;Lee and Johnson 2005).Light was supplied by

Table 3 .
Comparison of the functional response models of Neogobius melanostomus and Proterorhinus semilunaris at two temperatures; dAICc = difference in AICc value from the most parsimonious model (in bold); df = degrees of freedom; weight = Akaike weight.

Table 4 .
F and p values of the analysis of deviance for effects of temperature and predator species on maximum feeding rate.Bold values represent significant explanatory variables.