Modeling species invasions using thermal and trophic niche dynamics under climate change

Introduction

Global ocean temperatures have risen by 0.11°C per decade over the period 1971–2010 (Jones et al., 2013). This tendency is predicted to continue in the future, in particular for the Mediterranean Sea where warming scenarios project a mean annual temperature rise 1.25 times higher than the global average long term predictions for 2081–2100 (Jones et al., 2013). Changing marine temperatures play a role in modifying the distributional ranges and density of natural populations. The poleward shift of thermal habitats makes previously inhospitable areas accessible to more thermophilic (warmer-water) species, and makes the same areas less suitable for the native ones (Walther et al., 2002, 2009; Parmesan, 2006; Bazairi et al., 2010). These two processes are expected to open the so-called “invasion window,” changing an ecosystem's susceptibility to invasion (Carlton, 1996; Drake et al., 2006; Caplat et al., 2009). Climatic changes thus ultimately influence biodiversity (Bianchi and Morri, 2000), increase the risk of extinction for less adaptable species (Bazairi et al., 2010) and affect fisheries (Cheung et al., 2010). An increasing trend of non-indigenous thermophilic species has already been recorded in the Mediterranean Sea (Ben Rais Lasram and Mouillot, 2009; Zenetos et al., 2012), possibly resulting in an irreversible marine ecosystem shift (Ruiz et al., 1997).

Thermal and trophic opportunities, disadvantages and competition within ecological niches all play a role in determining a species' potential for successful invasion. Successfully invading thermophilic species are expected to take advantage of increased average temperatures (Dulčić and Grbec, 2000; UNEP - MAP - RAC/SPA, 2008; Dulčić et al., 2010, 2011; Brotz and Pauly, 2012; Zenetos et al., 2012; Pecarevič et al, 2013) but they also need to either find or make trophic niche space, such as by outcompeting native species (e.g., Pranovi et al., 2003). Meanwhile, some local native species might be negatively affected by higher temperatures, but may be able to exploit new trophic opportunities as a trade-off. It is therefore not trivial to predict the potential invasion of non-indigenous species (NIS) and their impacts (Jones et al., 2013; Pinnegar et al., 2014).

In this work we explore which are the potential thermophilic invaders and their ecosystem effects by representing combined thermal and trophic niches for both native and non-indigenous species. Dynamic thermal and trophic niches have been described using a trophodynamics model wherein water temperature is simulated to rise according to global and local sea surface temperature (SST) scenarios (Somot et al., 2006; Jones et al., 2013).

We used the northern Adriatic Sea (NAS) basin as the ideal hotspot to investigate the effects of climate change. This semi-enclosed and shallow basin (Figure 1) is exposed to cold north-easterly winds and receives cold waters from many alpine rivers that contribute to create a sub-atlantic climate, and therefore hosts several species adapted to boreal conditions (Tortonese, 1964). Moreover, its position within the Mediterranean, makes the NAS a cul-de-sac for these species, preventing further poleward migration as a response to temperature rise (Ben Rais Lasram et al., 2010). Multiple gears target several marine demersal and pelagic species living in the shallow and highly productive grounds of NAS (Pranovi et al., 2001; Pranovi and Link, 2009) with documented historical changes and effects on the marine communities (Fortibuoni et al., 2010; Barausse et al., 2011).

Figure 1. Map of the Northern Adriatic Sea. (A) Localization of the basin within the Mediterranean Sea; panel (B): detail of the area highlighting the main features determining a subatlantic climate, i.e., the north-easterly cold winds and the main rivers discharging in the Northern Adriatic Sea the cold waters from northern alpine areas.

Within this context, some questions arise in relation to the potential arrival of thermophilic non-indigenous species (NIS), such as:

– Are changes in community structure simply linear with temperature, implying easier predictability? How uncertain are our estimated biomass changes due to temperature changes?

– What is the combined effect of trophic and thermal opportunities in the dynamically changing environment both for native and non-indigenous species? Can we characterize species that benefit from climatic changes as those that show biomass increases (winners) and species harmed by climatic changes as those that show biomass decreases (losers) for both native and non-indigenous species?

– What role does the adaptation of native species to rising temperatures, changes in primary production, and changes in fisheries pressure play on the changes in community structure?

– Overall, how does invasion by NIS restructure the biological community?

We explore emerging ecosystem changes using a dynamic food web model for the NAS that:

(i) Includes, as initial conditions “seeds” for NIS at extremely low biomass;

(ii) has defined thermal niche preferences and trophic interactions for non-indigenous and native species, or groups of species;

(iii) is subjected to scenarios of water temperature increase.

With such a model structure we intend to simulate the potential opening of the invasion window due to climatic changes. The general goal is to explore the possible broad community consequences of the establishment of new species under climate change.

Materials and Methods

The Food-web Model

An existing model representing the NAS ecosystem in the 1990s (Zucchetta, 2002; Pranovi and Link, 2009) was updated with 2007–2008 data. The ecosystem was modeled with a mid-complexity structure using 30 living functional groups and 2 non-living groups (“Detritus” and “Carcass”), in order not to compromise the capability to represent the main trophodynamic processes, while maintaining an acceptable level of synthesis. Functional groups were defined by aggregating species with ecological similarities (Christian and Luczkovich, 1999) such as food preferences, predators, life history traits and habitat preferences (Table 1 and Supplementary Materials, Table S1). Some taxa were described by dedicated groups because of their importance as target species for local fisheries, namely Thunnus thynnus (Bluefin Tuna, BFT group), Platichthys flesus, Scophthalmus maximus, Scophthalmus rhombus, and Solea solea as a Flatfish group (FFS), Aequipecten opercularis and Pecten jacobeus as a Pectinidae group (PEC), and the venus clams Chamelea gallina and Venus verrucosa as a Veneridae group (VEN) (Table 1). Fishing activity was represented through five fleets: tuna fisheries (that include longline and purse seine), mid-water trawl, otter-trawl, hydraulic dredge, and rapido trawl (a typical fishing gear used in the Adriatic to harvest flatfish and scallops—see Pranovi et al., 2001). These fisheries are modeled by representing both commercial catches (landings) and by-catch of non-commercial species: this latter constitutes the discard, which makes up a particularly large fraction of catches for otter and rapido trawl (Pranovi et al., 2001). Therefore, the abiotic compartment “Carcass” has been added to account for the discard from fishing activities, and the consequent scavenging processes that are important in the NAS ecosystem dynamics (Zucchetta, 2002; Pranovi and Link, 2009).

Table 1. Functional groups included in the model, with input and output (in bold) basic parameters of the Ecopath mass balance: trophic level (TL), production rate (P/B), consumption rate (Q/B), Ecotrophic efficiency (EE), growth efficiency (P/Q).

Dedicated functional groups, exclusively composed of thermophilic non-indigenous species, were introduced at all trophic levels as potential homologous of some native ones. The invasive groups represent species recorded recently in the basin such as the planktivorous fish Sardinella aurita and the barracuda Sphyraena viridensis, or NIS recorded in other parts of the Mediterranean and whose arrival in the NAS is likely to occur in the near future (e.g., the opistobranch Halgerda willeyi, Zenetos et al., 2012). Parameters for the new groups were calculated taking into account productivity and consumption rates of the species composing the group, as well as the specific diets, which were derived from the data available on SeaLifeBase (Palomares and Pauly, 2014). The initial biomass for the invasive groups was set to very low values, so to represent a “pre-invasion” status.

The modeled functional groups and their main parameters are summarized in Table 1, while Table S1 in the Supplementary Materials details the sources used to determine the parameters, and Table S2 shows the diet composition matrix. A synthetic view of the model structure and main energy flows are reported in the Figure 2.

Figure 2. Scheme of the structure of the food-web model. Functional groups and fisheries are arranged clockwise in decreasing order of trophic level (TL). Barplots of the outer circle represent optimal temperature and maximum temperature (T_opt and T_max values in cyan and red bars, respectively) in the scale 0–30°C for the functional groups, which are represented in the inner circle as labeled boxes (group codes are reported in Table 1). The width of the boxes is proportional to the group biomass, and the color indicates the group type (Green = Native, Orange = Invasive, Blue = Fishery, Gray = Non-living). Links, representing trophic interactions, are proportional to the flow between two groups; inbound flows are those closest to each box.

The model, built in Ecopath with Ecosim 6.4 (Christensen et al., 2000; Christensen and Walters, 2004; Christensen, 2009), is based on a system of ordinary differential equations, describing the biomass variation of each functional group over time (t) as follows:

where B_i is the total biomass of the ith group composing the modeled food-web, (P/Q)_i is its production/consumption ratio, Q_ji is the consumption of group i (predator) on group(s) j, Q_ij is the consumption by group j on group i (prey), I_i is the immigration rate, M_i is the natural mortality rate, F_i is the fishing mortality rate and e_i is the density dependent emigration rate. Initial conditions and parametrization were set on the basis of the mass balance Ecopath model (Table 1). In the NAS model e_i and I_i were set as zero, M_i was estimated from Ecopath production rate (P/B)_i and ecotrophic efficiency (EE_i) as M_i = (1-EE_i)(P/B)_i. Initial estimates of fishing mortality were calculated from the ratio between the catches [Y^E_i = Y_i(t = 0)] and biomass [B^E_i = B_i(t = 0)] in Ecopath, i.e., F^E_i = Y^E_i/B^E_i. A time series multiplier ff(t) was entered in Ecosim to modulate changes of fishing mortality over time, F_i(t) = ff(t)* F^E_i. It has to be noted that in EwE, non-living compartments (Detritus and Carcass) are set to be flow-based so their dynamic is almost insensitive to their initial biomass.

Each element Q in Equation 1 is based on the foraging arena (Walters et al., 2000; Walters and Christensen, 2007; Christensen, 2009; Ahrens et al., 2012), which assumes that only a portion of the prey can be vulnerable to predator: this accessible biomass is regulated by a vulnerability rate (v_ij) and the quantity of prey consumed by a predator will be:

where a_ij is the search rate and the forcing function f(T,t) can be used to account for external drivers changing over time, such as temperature T (Ainsworth et al., 2011). Ecopath mass balance biomasses B_i and B_j, predator consumption rate (Q/B)_j and diet proportion DC_ij, are used to estimate a_ij. More details on Ecopath and Ecosim parametrization can be found elsewhere (Christensen and Walters, 2004b).

For consumer groups, to relate the growth rate of a population to temperature, we set the forcing function as the one initially proposed by (Lassiter and Kearns, 1974):

where T is the SST, T_opt is the optimum temperature for a given species, at which, with other factors constant, the consumption and growth is maximum; T_max is the maximum tolerated temperature, above which the consumption and growth is null, and c is a parameter accounting for the sensitivity of a species' consumption (i.e., its productivity) to temperature variations.

For primary producers, two types of forcing functions were considered to modify productivity, a linear variation with temperature, and an exponential one (Eppley, 1972):

Estimating Temperature Dependence Parameters

Although, some empirical data regarding thermal preferences and effects of temperature on physiological processes are available in literature or datasets (e.g., Coutant, 1977; Block et al., 2001; Roessig and Woodley, 2004) few of them regard Mediterranean marine species and potential invaders. Therefore, in order to be consistent across the many species that we included in the model, temperature dependence parameters were estimated using a geographical distribution approach similar to what has been adopted in other cases (Dulvy et al., 2008; Cheung et al., 2013; Parravicini et al., 2015).

The parameters T_opt, T_max, and c were estimated for each species on the basis of their latitudinal distribution in the northern hemisphere and by reconstructing a mean latitudinal SST profile to relate latitude to temperature. Species latitudinal distribution was derived from Global Biodiversity Information Facility database (www.gbif.org) considering occurrences up to year 2000 to exclude possible recent range modifications due to climate change. A yearly mean latitudinal SST profile for the northern Atlantic Ocean was reconstructed using SST data from 1945 to 2000 of the online NOAA database (nodc.noaa.gov). The T_opt for each species was defined as the mean annual SST at the latitude corresponding to the median of the species' latitudinal range. Similarly, T_max was defined as the mean annual SST at the latitude corresponding to the 5th percentile of the species latitudinal distribution. The parameter c, that indicates the rate of change of a group's consumption and productivity with changes in SST, was applied considering a species preferred habitat and its possible relation to the surface. Thus, high values of c were assigned to pelagic species, intermediate values to demersal ones, and lowest values of c were assigned to benthic species in order to represent decreasing response of species to SST changes. This parametrization permits to represent the fact that a species that can move and adjust its location quickly will maintain the highest productivity relative to its specified base rate.

Once estimated for all species, the temperature parameters for each functional group were calculated using biomass weighted averages of the values of the species composing the group. Our resulting temperature parameters were compared also with values of T_opt and T_max from (Cheung et al., 2013) for validation (see Table S2, Supplementary Material).

Sensitivity Tests and Scenario Analysis

The model was tested for sensitivity to SST trend changes, and to the parameters of temperature forcing function in Equation 3.

The initial temperature was set to the mean SST registered in the basin in 2007 (T_start = 18.4°C); the final temperature change was comprised between +0.8 and +1.2°C in 10 years, and then kept constant for 20 years to reach steady state (T_end = 19.2 ~ 19.6°C). A set of 20 simulations over 30 years was run each with a different T_end at regular intervals of 0.02°C. For each simulation, the final biomass (B_end), the relative biomass variation (B_end/B_start) and the value of the forcing function (f(T)) were calculated for each group; further, the relative sensitivity, S (Campolongo et al., 2000), was calculated as follows:

where ΔB = B_end – B_start and ΔT = T_end – T_start. This relative sensitivity is a measure of expected changes around the starting conditions (local sensitivity) and it indicates the expected relative change in biomass (% change) due to a unit relative change in temperature (i.e., 1%) (Zádor et al., 2006).

Moreover, in order to test the sensitivity of the model to the uncertainty on temperature parameters, 100 simulations were performed assuming a linear temperature rise of +1°C in 10 years and using different combinations of T_opt, T_max, and c for every functional group. Parameter values were randomly extracted from normal distributions centered on the actual estimated parameter and having SD = 0.7. The distribution of results in terms of relative biomass variation (B_end/B_start) was used to evaluate the uncertainty of predicted changes.

In order to test ecological changes induced by adaptation, changes in primary productivity (PP) and fisheries and their combined effects a set of 24 scenarios (Table 2) grouped into six groups were implemented. The reference scenario with current PP and fisheries pressure was run for 50 years under the IPCC A2 that implies a 0.03°C SST rise per year (applied for the first 30 years) on the basis of temperature projections by Somot et al. (2006).

Table 2. Conditions and drivers used in the scenario analysis on top of the reference climatic change.

Then the capability of autochthonous groups to adapt to temperature changes was represented by applying the temperature-based forcing functions to the search rates of every group in the system (non-adaptation; scenarios 1–3) or to all groups but with a delay for native groups (adaptation, scenarios 4–6) on the assumption that the search rate of the native is affected by T rise only after a few years after T starts changing (Table 2). Specifically, a delay of 5 years was assumed for groups with Topt <= 10°C (FFS, PEC, VEN, MHR), and a delay of 10 years was applied to all the others.

PP was represented by (i) no change (scenarios 1 and 4), (ii) an exponential increase in production with temperature (scenarios 2 and 5), and (ii) a linear decrease of 20% in 30 years representing the oligotrophication trend actually detected in the basin (Giani et al., 2012; scenarios 3 and 6) (Table 2).

Moreover, the reference simulation (current fishing effort exerted on native species only; scenarios C) was compared with simulations representing as effort reduction to zero (no fishing; scenarios A), an effort reduction by 20% (fishery on all groups; scenarios B), and the new invasive species as an opportunity for fishermen (F for invasive species equal to the native counterpart; scenarios D) (Table 2).

Ecosystem Structure, Functioning, and Ecological Niches

Ecosystem effects of invasion in combination to other factors represented under different scenarios were evaluated by looking at changes in biomass of functional groups and by analysing the Kempton's Q₇₅, an index that synthesizes the community structure in the form of biodiversity index (Ainsworth and Pitcher, 2006). The Q₇₅ index is calculated as:

where S is the total number of functional groups and R₇₅ and R₂₅ are the biomasses at the 75th and 25th percentile, respectively. Moreover, for distinguishing the prevalence of trophic effects or of temperature-related effects on the various groups, the Mixed Trophic Impact (MTI) index was used (Ulanowicz and Puccia, 1990). MTI results from the application of the input-output analysis on the food web flows and quantifies the direct and indirect trophic effects that one group has on each other of the food web (See Supplementary Materials, Figure S1). The sum of MTI by row represents the overall effect (MTI_j) of a group on the food web and it is a prediction of biomass changes as a consequence of the biomass variation of all the other groups in the system (Libralato et al., 2006). Therefore, by plotting biomass changes obtained from simulations of increasing SST against MTI_j (that do not include temperature effects) enables to distinguish results that are fairly in line with MTI_j predictions and that are therefore mainly driven by trophic (direct and indirect) opportunities and disadvantages, from results that contrast with MTI_j, for which direct and indirect effect of temperature rise can be considered to be dominant.

Results in terms of biomass change with SST change scenarios allows us to distinguish between groups with positive (winners) and negative (losers) biomass change: these were compared with predator and prey niche overlap (Christensen et al., 2000) to identify possible trophic niche shifts and partial niche substitutions.

Results

The analysis of sensitivity to the variability of T highlights the existence of nonlinear effects within the modeled food web (some example functional groups are represented in Figure 3). In particular, analyzing the relationship between the value of f(T) (Equation 3) calculated for T_end and the biomass (Figure 3A) or the relative sensitivity index S (Equation 5, Figure 3B), nonlinear changing rates emerged for every functional group.

Figure 3. Relationship between T increase and relative biomass change (A), Temperature forcing function (see Equation 3) and relative biomass change (B), T increase and relative sensitivity as from equation 5 (C) and temperature forcing and relative sensitivity (D). For the sake of clarity, results are reported only for some representative groups, i.e., the BFD, Benthic Feeders group; MMF, Macrobenthic Mixed Feeders; NFI, Nekton Feeders Invasive; MFI, Macrobenthic Filter Feeders Invasive.

The results of uncertainty tests on the parameters of temperature equation by species (Equation 3) showed considerable changes and variability in the functional group biomass estimated by the model after 30 years of simulation (Figure 4). Although, as expected, all invasive species benefited of environmental changes, the resulting positive biomass change was different for each functional group. Moreover, uncertainty analysis showed that new ecosystem conditions might have non-zero probabilities to result in negative effects for some invasive thermophilic groups such as MOI, MMI, JEI, and MZI (Figure 4). Some native groups showed large probability of increasing biomasses (BFT, PLT, MEI, JEL, MIZ, BPL), whereas others (NFD, CPH, FFS, MOP, PEC, VEN, MZP, MHR) had higher probabilities of decreasing their biomass in the new environmental conditions. Some groups, such as BFD, MMF, MFF, and MDT, had very large probabilities of maintaining current densities even with new climatic conditions. Variability rarely exceeded ±30% of the median value (Figure 4). In general, medium-high TL groups showed larger dispersion of results than low TL groups. Moreover, although MOI and MMI had the largest variability in final biomass, several native groups (BFT, BFD, PEC, VEN, MHR) were more sensitive to simulated changes than any other invasive groups.

Figure 4. Sensitivity tests on the temperature forcing function parameters [f(T,t), Equation 3]. Results are reported in terms of biomass at the end of the simulations (Bend) relative to initial biomass (Bstart). Box plots represent the median (black line) and interquartile range (gray box) for relative biomass of each functional group. Native groups are reported on the left, invasive groups on the right. The solid line (Bend/Bstart = 1) represents no change. See Table 1 for explanation of functional groups codes.

Comparing biomass changes for both native and invasive functional groups for simulations based on the SST increase scenarios with simulations where the SST increase is combined with modifications of the primary production and fishery (Figure 5) permits us to disentangle climate from the other sources of perturbation. Results seem to confirm the heterogeneity in terms of responses of different groups to the SST increase, with positive, negative and neutral “behavior” (Figure 5, top left panel). Given a fishery scenario, the increase or decrease in PP produced a unidirectional shift of all groups toward higher or lower biomass (Figure 5, panels in central and right columns respectively) compared to constant PP (panels in the left column). For the current fishing conditions (lowest row of panels, D), this translates in a shift above and below the reference 1:1 line. Fishing seems to produce a strong effect on the ecosystem: current fishing conditions including (row D) or excluding (row C) fishing on invasive species have prevailing and dominating effect on biomass changes, regardless of PP changes. Reducing (row B) or eliminating fishing (rows A) implies considerable restructure of the food web due to the biomass increases of both invasive and native functional groups.

Figure 5. Results of the simulations based on simultaneous SST increase, changes in primary production and fishing activities. Results for all scenarios are reported in terms of relative biomass variation with respect to initial conditions for each functional group. In all these scenarios, the native species are allowed to adapt to the SST increase, at least in the first phase of the simulation (analogous scenarios with no adaptation are reported in the Supplementary Material). According to the scenario definitions in Table 2, A = no fishery, B = fishery on all groups, but reduced to 20% of the present effort level, C = fishery only on native species, D = fishery on all species at the present effort level. Functional group codes are defined in Table 1.

Changes in the ecosystem structure, due to thermophilic species invasion, were well captured by the Kempton's Q₇₅ index, showing a common pattern that supports the ecosystem development theory (Odum, 1969): the initial steady rise due both to the increase of invasive biomass and initial reduction of the most abundant native species, is followed by a maximum plateau (generally, attained before the forcing function stops acting, vertical dotted line) and then a drop, after which the index tends to stabilize at a value higher than the initial one (Figure 6). It is worthy to note that the adaptation of native species (scenarios 4–6) produced a more gradual variation, with a reduction and shift of the peak. Changes in PP slightly affected the peak (both in term of position and value; scenarios 2, 3, 5, and 6). Finally, modifications in the fishing pressure regime exacerbated the effects in terms of changes in the community structure, amplifying the peak and shifting it in the time and leading to ultimate reductions in biodiversity (Figures 6C,D).

Figure 6. The Kempton's Q75 index, estimated for different tested scenarios. Reference scenario with current primary productivity PP and fisheries pressure was run for 50 years under the IPCC A2 that implies a 0.03°C SST rise per year (applied for the first 30 years) on the basis of temperature projections by Somot et al. (2006). Scenarios 1–3 are simulated with no adaptation of native species to changes in SST, while in scenarios 4–6 adaptation to SST is simulated. Scenarios 1 and 4, represent no changes in PP, while scenario 2 and 5 increase, and 3 and 6 decrease of PP. (A–D) refers to scenarios with increasing fishing pressure: (A) no fishing; (B) actual effort reduced by 20% in 5 years; (C) only resident species are targeted at actual effort; (D) all species are targeted at actual fishing effort. All scenarios are summarized in Table 2. The vertical dashed line represents the point beyond which the forcing stops acting.

By comparing the MTI with the biomass changes resulting from the sensitivity tests, it was possible to highlight the model capabilities to distinguish the prevalence of climate-related or trophic-web related effects. Plotting the relative biomass variation index against the MTI index (Figure 7), different groups can be detected:

• (A) those for which trophic effects were indicated to be negligible (MTI < 2), but the high biomass variation [positive or negative, (B_end/B_start) > 1.5 or <0.5], highlighted that their changes were driven mainly by temperature (e.g., FFS, MHR, MHI, PLI, MDI);

• (B) those with little/intermediate biomass variation 0.5 < (B_end/B_start) < 1.5, but with large expected trophic effects, e.g., MTI >> 2 (BFT, MMI, CPH), for which it seems that temperature and trophic effects act in an antagonistic way;

• (C) groups with high biomass change [(B_end/B_start) > 1.5 or < 0.5] and high trophic effects (MTI >> 2) were indicated (NFI, MOI, PEC), for which a synergistic effect of trophic and temperature impacts can be hypothesized;

• (D) those groups with low intermediate changes in biomass and low expected trophic effects that showed mixed somewhat small effects (most functional groups).

Figure 7. Relationship between relative biomass variation (B_end/B_start) from reference simulation (whiskers indicate range of variability from sensitivity tests) and overall impact index (MTI_j) for all the trophic groups. Native groups are in blue and invasive groups are in red. The solid line (B_end/B_start = 1) defines winners (above the line) and losers (below the line) in a context of climate change. Dotted lines separate areas of the graph with different dominance of trophic and thermic effects.

This analysis (Figure 7) allows discriminating between trophic and thermal drivers and resulting identification of winners (biomass change >1) and losers (<1). Combining this analysis with the niche overlapping index permits us to describe possible trophic substitutions and adjustments in the biological community (Figure 8). High values of index of predation overlap between losers and winners (dark blue cells in Figure 8) suggest, for example, that NFI and PLI (invasive winners) can substitute as predator NFD, FFS, and CPH (native losers); the invasive MOI potentially substitutes (at least in part) the native groups MOP and MDT as predator. As prey, the invasive NFI substitutes the native NFD, invasive PLI takes advantage of decreasing VEN and PEC, and MOI substitutes FFS (Figure 8).

Figure 8. Prey and predator niche overlap (lower and upper diagonal part of the matrix, respectively). Intensity of the color represents strength of the niche overlap (light = low niche overlap; dark = high niche overlap). Groups (codes are explained in Table 1; invasive groups are bold) are ordered in terms of relative biomass change (y-axis of Figure 7), thus to distinguish losers (A) from winners (B, see also colors). The high value of niche overlap for winners and losers permits the identification of niche substitution. For example, NFI is a winner and as prey, substitutes mainly the niche of the loser NFD (lower diagonal matrix). As predator, instead (upper diagonal matrix) the winner NFI occupies the niche of a few losers such as FFS, CPH, and NFD.

Discussion

Changes in species distributions is a climate change effect of primary interest and can integrate changes in phenology, condition, and biological communities (e.g., Parmesan, 2006; Lejeusne et al., 2010). However, previous analysis already showed that thermal niche alone cannot completely explain invasive species distributions (e.g., Parravicini et al., 2015). Combining thermal and trophic niche, for both native and non-indigenous species, in the food web model resulted a successful approach to represent the invasion and explore ecological effects of NIS. Niche availability, indeed, has been demonstrated as one of the most important aspects in determining the invasion success of NIS (Ward and Master, 2007). The so-called “invasion window” (Carlton, 1996; Drake et al., 2006; Caplat et al., 2009), depends on the niche availability, i.e., the set of ecological conditions that can transform an invasive species into “the right species at the right moment” (Pranovi et al., 2006).

Climatic changes, and in particular ocean warming, are expected to act as a factor contributing to the opening of the invasion window beyond the direct favoring of thermophilic species. In fact, thermal regime modifications that negatively affect native species might destabilize community structure, increasing niche availability and so further favoring the success of invasive species establishment (Wolkovich and Cleland, 2014). In this context, an explicit and “open” modeling scheme was developed instead of the ordinary approach of comparing characteristics of marine food web models with and without the addition of one or more invasive species (Arias-González et al., 2011; Pinnegar et al., 2014). Here, the SST was used as a forcing function driving the ecological dynamics of the different groups in the model, which includes a “seed bank” of non-indigenous species, thus permitting us to represent the invasion window opening and the establishment of invaders to emerge as a combined result of both direct and indirect trophic and thermic opportunities. Results of simulations representing an increase in the range of 0.8–1.2°C over 30 years and accounting for the uncertainty in the biological thermal preferences, robustly highlighted complex non-linear and non-homogeneous changes across groups of species. The establishment of invaders according to better thermal and trophic fitness, in fact, was accompanied by highly variable effects (biomass changes) on native species. Overall the biomasses of some species with boreal affinity, i.e., more sensitive to increasing temperature, such as the flatfishes (FFS), were importantly reduced, and the trophic overlap index permitted us to see that such trophic niche depletion provides opportunities for invasive groups even if taxonomically and ecologically very distant, such as the Planktivorous fish Invasive group (PLI, Figure 8). This exemplifies the complex processes involved in the reorganization of the entire food web due to community structure modifications that can be manifested by combining trophic with thermal niche changes. This increases the number of dimensions explicitly accounted for, thus approximating the real world multidimensional hypervolume of niche space (Pianka, 1981; Rosenfeld, 2002; Okey, 2003).

The combination of increase/decrease of species in relation to their thermal affinity is expected to produce a more or less gradual (depending on the rate of warming) substitution of temperate and boreal species with warmer-water species, as constrained also by variables such as pH, dissolved oxygen, and other habitat variables. Poleward shifts in distributional areas were predicted (Albouy et al., 2012) and already recorded (Parmesan, 2006; Cheung et al., 2013), all around the world. According to the shifting scheme just reported, every species can simultaneously act as native loser and non-indigenous winner, depending on the perspective (original or new ecosystem, respectively), and the phenomenon of native species substitution with thermophilic non-indigenous invaders (eulerian view) can be seen as a simple species shifts (lagrangian view point). However, at the “end of the line,” i.e., where no further Northward shift (in the boreal hemisphere) is possible due to the cul-de-sac effect, or other expressions of poleward boundedness, the destabilization and resistance of the native community that cannot shift may be more evident, leading to a more complex opening of invasion windows.

This is the case of semi-enclosed basins, such as the Mediterranean Sea and its northern sub-basins, or the high latitude systems, as for example the Barents Sea, which act as a cul-de-sac for the native species (Ben Rais Lasram et al., 2010), and prevent further migration. The direct and indirect effects of warming and invasion estimated using the food web model for the Northern Adriatic Sea, never resulted in a complete overturn of a group of native species or in a bloom of an invasive species, at least within the projected perturbations. Possibly the cul-de-sac regions represent areas where climate-induced changes are more complex than a simple poleward shift. These areas, therefore, provide good case studies to explore the ecological dynamics related to the simultaneous change of the environmental conditions (due to water warming) and of community structure (due to species invasion and native species contraction).

Invasion by thermophilic NIS implies important changes in the marine ecosystem structure that are confirmed by our results in terms of changes in biomass of ecosystem components after the establishment of non-indigenous species. However, the results obtained in the present study showed that the direction of change is not always as expected, as indicated by the increase of the diversity index (after a transient period with large increase of diversity).

The consistent increase of diversity, both under changes in productivity of the system and adaptation of the residents, was not only due to the increase in biomass of invasive groups but primarily by the reduction of biomass of residents, which increased the evenness among groups. Despite a complex set of biomass changes, the fact that the diversity index stabilizes at increased values after invasion indicates that in such enclosed systems, the establishment of invasive species represents an enrichment in ecosystem structure. This is not necessarily a positive change since the increase of biodiversity might be accompanied by loss of non-trophic function or habitat (Galil, 2007) that are not included in the analysis here. Overall the absence of complete species substitution suggests a possible contribution of resident species in increasing the resilience of the system to invasion. As highlighted by Carey and Wahl (2010), indeed, the biodiversity of native species can mitigate the impact of NIS, for instance by saturating niche availability.

Within this context, results obtained at different fishing pressure illustrate how much fisheries can destabilize ecosystem structure and processes in terms of dynamics of the diversity index. In particular current fishing pressure (only on native species; or on both native and invasive species) represents a constraint stronger than invasion for density changes and it has a sort of “homogenization effect” on almost all the groups, strongly diminishing effects of both SST and primary production changes. In contrast, high fishing pressure on native or both on native and invasive influenced the synthetic diversity index by simply delaying the effects of invasion (Figure 6). Reducing fishing pressure in the context of invasion, instead, resulted in a complete restructure of the food web, with significant biomass changes (Figure 5) and with limited effects on the diversity index (Figure 6).

Overall these results indicate that invasion affects both the structure, with possible effects on the functioning of the ecosystem. Adaptation of native species has minor effects on density but results in delaying dynamics of the diversity index, while fisheries substantially impact species biomass changes and thus ecosystem structure.

The combination of thermal and trophic niches in the food web model of species invasion suggests that, due to direct and indirect effects, it is fairly difficult to predict winners and losers a priori. Although for some species the effect of temperature is dominant (such as FFS), while others have a dominant trophic effect (such as BFT), the majority of modeled groups indicate that they are affected by a combinations of the two factors, with varying interactions of the two effects. This is confirmed by the fact that in only a few cases the winners occupy similar niche of a loser.

In conclusion, the results indicated the utility of an enhanced and explicit ecological niche approach for advancing in the analysis of possible effects of species invasion. Moreover, results highlighted the complexity of dynamics linked to an invasion by thermophilic non-indigenous species, not only in terms of predicted impacts and their strengths, but also in terms of the directions of resulting changes in biomass.

Within the context of habitat modifications induced by climate change (e.g., SST increase), the thermophilic NIS invasion is generally viewed as a gradual substitution between native and invading species, both shifting northward (Cheung et al., 2010). In semi-enclosed basins, however, the two groups are forced to compete, with the native species able to partially counteract (resist) the invasion, or at least persist after invasion, as indicated also by the predicted delay of native species adaptation to the SST increase.

These findings highlight the importance of native species for ecosystem resilience, facilitating the system's ability to cope with the modifications induced by climate changes. In this respect, fishing operations play an important role in shaping community structure as it was also emerging from previous analyses of historical data (Fortibuoni et al., 2010; Barausse et al., 2011). In light of the fact that anthropogenic pressure (here exemplified by fisheries) is the main stressor responsible for the loss of ecosystem resilience (Fortibuoni et al., 2015) to climate-related changes, the process of defining management strategy that policy makers and managers would normally consider are those that reduce fishing pressures that disproportionately degrade the system, or reduce the resilience. However, the present approach could be useful for evaluating management strategies that include deliberate and profitable counteraction of climate effects, such as targeting newly arriving species in a way that maximize or preserve the historical values or other values in the overall social-ecological system.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors wish to thank Elisabetta Betulla Morello (Marie Curie fellow, CNR – ISMAR, Ancona, Italy), Tomaso Fortibuoni (ISPRA, Chioggia, Italy) and Thomas A. Okey (School of Environmental Studies, University of Victoria, Canada) for providing valuable comments and suggestions on an early draft of this paper. SL has been partially funded by the Flagship Project RITMARE - The Italian Research for the Sea - coordinated by the Italian National Research Council and funded by the Italian Ministry of Education, University and Research within the National Research Program 2011-2013.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fmars.2015.00029/abstract

References