Highlights

  • Degradation on HFS induced by OA has unexpected cascading effects on marine ecosystem functioning and services.

  • Food web model allows the quantification of a decrease in the system size and a reorganization of energy flows in the coastal ecosystem due to OA.

  • Scenarios with the integrated modeling approach show that all the proxies for ecosystem services significant decrease in their values.

Introduction

Global changes are resulting in progressive modification of relevant oceanographic characteristics, such as ocean temperature and currents, and the acidification of global oceans (IPCC 2019). How marine ecosystems may respond to such environmental changes and how ecosystem functioning will be impacted are far from completely understood. Given the importance of goods and services that marine ecosystems provide to humankind, it is relevant to evaluate the implications of global modifications in ecosystem functioning for socioeconomic systems and human well-being (IPBES 2019).

Although the comprehensive evaluation of all the effects induced by climate change is an overwhelming task, the quantification of effects induced by a specific pressure on key sensitive species and related cascading effects might provide an opportunity to quantify the implications of such relevant global impacts. Here, we aim to explore and quantify the potential effects of habitat-forming species HFS loss due to OA on marine population, community and ecosystem using a modeling approach.

Our study focuses on the Mediterranean Sea, which is a basin that is considered very sensitive to global changes (Giorgi 2006) and to a variety of intense human pressures (Micheli and others 2013) and might be considered a miniature ocean, since most of the processes that occur in the global ocean also take place here but at smaller and faster space/time scales (Bethoux and others 1999; Crise and others 1999; Lejeusne and others 2010).

We focus on two Mediterranean marine habitats that are important suppliers of highly valuable ecosystem services (ESs), play a fundamental role in supporting human well-being and are considered to be vulnerable to anthropogenic pressures and ocean acidification (OA): coralligenous reefs and Posidonia oceanica meadows (see the review in Zunino and others 2019).

These are considered to be habitat-forming species (HFS) that are able to create complex environments sustaining the functioning of ecosystems and complex food webs (D’Onghia and others 2010), representing hotspots of biodiversity providing refuge and nursery habitats for many commercially important species. For example, coralligenous algae and corals are preferred habitats for the dusky grouper, Epinephelus marginatus, which is classified by the IUCN as “endangered” (Cornish and Harmelin-Vivien 2006). Changes in habitat complexity due to a decrease in HFSs represent a key type of alteration for benthic systems, with potentially large effects on ecosystem functions and, consequently, on a cascade of ecosystem services (Sunday and others 2016).

HFS degradation and biomass declines were described in a recent analysis that qualitatively defines the ecosystem impacts of the loss of P. oceanica meadows and coralligenous habitats due to the degradation effects of OA (Zunino and others 2019). Qualitative impacts are quantitatively parameterized and translated into food web models.

Ecosystem models play an instrumental role in quantifying ecological relationships in the environment (Rose and others 2010; Sarmiento and Gruber 2013), permit the analysis of processes involving different ecosystem components in a single coherent framework (Sala 2004; Mackinson and others 2009; Libralato and Solidoro 2009). Ecosystem models allow the assessment of the cumulative effects of direct and indirect interactions among multiple species (Peck and others 2016; Celić and others 2018) and might result in tools for testing management strategies (Essington and others 2017).

The analysis is conducted by extending an ecosystem model representing a well-studied coastal marine protected area in the NW Mediterranean Sea (Prato and others 2016). The model allows the representation of scenarios of increasing levels of OA impacts on HFS, thus permitting evaluation of the ecological cascading effects and the long-term sensitivity of the ecosystem to OA. Furthermore, the changes in fishing captures and opportune indicators (for example, Plummer and others 2013) derived from the modeled scenarios allow a novel inference on the impacts of OA degradation on ecosystem services.

This work aims to model the impact of OA on HFS biomass and the cascading trophic and non-trophic effects (hereafter called “mediating effects”) on the biomass of ecosystem components, ecosystem processes, interactions, and properties. Moreover, this work quantifies the potential long-term effects and implications of OA for ecological, social, cultural and economic aspects related to marine systems using a set of proxies for ecosystem services. Thus, the study attempts to quantitatively assess the cumulative impact of OA at the whole-ecosystem level, considering ecosystem functioning and services. Despite the limited extent of the area represented, the analysis reveals patterns of change and ecosystem feedback that might have a more general validity.

Materials and Methods

The food web model (Prato and others 2016) was developed using the popular software Ecopath with Ecosim (Christensen and Walters 2004), and it represents the dynamic evolution of multiple interacting species in a marine protected area in the northwestern Mediterranean Sea under present conditions (years 2007–2014). The modeled area is mainly characterized by hard bottoms (51% rocky habitat, 31% coralligenous habitat), some P. oceanica meadows, and shallow sands (18% of the area overall).

The model includes 33 functional groups, where a functional group represents either a single species or a group of species (labeled with a “+,” in Table S1) linked by trophic relationships (Figure S1). The model includes three groups representing primary producers (that is, Posidonia oceanica, seaweed and phytoplankton); eight groups of invertebrates, including one cephalopod group and two for zooplankton groups; and detritus (Prato and others 2016). The model represents all the functional groups that are relevant for our purpose (that is, corals, P. oceanica meadows) and both as adults and juveniles of commercially important species (see also Prato and others 2016). The model also includes a representation of small-scale professional fisheries and recreational fisheries that are active in the area, represented with their portfolio of catch species. Data quality was assessed by Prato and others (2016) applying the food web diagnostics proposed by Link (2010) (see Figure S2).

The existing Ecopath model provided the initial conditions for the dynamic Ecosim routine, which describes the biomass variation of each functional group (Bi) over time (t) as follows:

$$\frac{{\text{d}B}_{i}}{\text{d}t}={g}_{i}\sum_{j}{Q}_{ji}-\sum_{j}{Q}_{ij}+{I}_{i}-\left({M}_{i}+{F}_{i}+{e}_{i}\right){B}_{i}\left(t\right)\;\text{for consumers}$$
(1a)
$$\frac{{\text{d}B}_{i}}{\text{d}t}={PP}_{i}\left(t\right)-\sum_{j}{Q}_{ij}+{I}_{i}-\left({M}_{i}+{F}_{i}+{e}_{i}\right){B}_{i}\left(t\right)\;\text{for producers}$$
(1b)

where dBi/dt represents the change in the biomass of group i occurring over a time interval of dt; g is growth efficiency; Qji defines the trophic relationships between a consumer (j) and a prey (i); Ii and ei are the immigration rate and emigration rate, respectively; Fi is the fishing mortality rate; and Mi summarizes all other mortality not explicitly represented in the model by predation and fishing mortality. PPi is the function that drives the primary production of autotrophs over time, as follows:

$${PP}_{i}\left({B}_{i}(t)\right)= \frac{{r}_{i}{B}_{i}(t)}{\left(1+{k}_{i}{B}_{i}(t)\right)} $$
(2)

where ri is the maximum production rate, and ki is a species-dependent parameter.

Food web dynamics are driven by trophic interactions based on the foraging arena theory (Eq. 3), in agreement with the rate of predation depending on the prey biomass Bi(t), the predator biomass Bj(t), the effective search rate for species i by species j, aij, and the rates at which species i moves between the vulnerable and invulnerable prey pools, vij (Christensen and Walters 2004):

$${ Q}_{ij}\left({B}_{i}\left(t\right),{B}_{j}\left(t\right)\right)= \frac{{a}_{ij}{v}_{ij}{B}_{i}{\left(t\right)B}_{j}\left(t\right)}{\left({v}_{ji}+{v}_{ij}+{a}_{ij}{B}_{j}\left(t\right)\right)} $$
(3)

Indeed, predator–prey interactions are assumed to take place primarily in restricted “foraging arenas,” where prey become vulnerable to predation (Walters and Christensen 2007; Ahrens and others 2012). The Ecosim module is run with the default vulnerability setting because of the lack of time series of data for fitting and tuning these parameters (Heymans and others 2007).

In this dynamic module, the OA impacts have been described by considering the following:

  1. (a)

    The direct change in a habitat-forming species, h, assuming that the degradation on HFSs can be represented as a reduction in their productivity, which results in an HFS biomass reduction, triggering cascading changes in the trophic interactions between species h, its predators and its prey, Figure 1a. Mathematically, this is parameterized by replacing the constant parameter ri in Eq. (2) with a function, ri S(pH(t)), which depends on the time-varying value of pH(t):

    $${PP}_{h}\left({B}_{h}\left(t\right),t\right)= \frac{{r}_{h}{ {S}_{h}\left(pH(t)\right)\cdot B}_{h}(t)}{\left(1+{k}_{h}{B}_{h}(t)\right)}$$
    (4)

    Similarly, for heterotrophs, the constant parameter aij in Eq. (3) is replaced with a function that depends on pH(t), ai,j*Sij(pH(t)):

    $${ Q}_{ij}\left({B}_{h}(t),{B}_{j}(t)\right)= \frac{{a}_{ij}{{S}_{ij}\left(pH\left(t\right)\right) v}_{ij}{B}_{h}(t){B}_{j}(t)}{({v}_{ji}+{v}_{ij}+{a}_{ij}{{S}_{ij}\left(pH\left(t\right)\right) B}_{j}(t))}$$
    (5)
  2. (b)

    The direct change (biomass reduction) in habitat-forming species h simulated as indicated above can trigger an additional change in trophic interactions that does not involve species h (see Figure 1b) but occurs among other species. Mathematically, this is parameterized by replacing the constant parameter aij, describing the selected trophic relationships between species i and j (different from h), in Eq. (3) with the function aij·Mij(Bh(pH(t)), which depends on the biomass of the habitat-forming species h. Analogously, the constant parameter vij is replaced by the function vij·Mij(Bh(pH(t)):

    $${Q}_{ij}\left({B}_{i}(t),{B}_{j}(t)\right)= \frac{{a}_{ij}{{M}_{ij}\left({B}_{H}( pH(t))\right)v}_{ij}{B}_{i}(t){B}_{j}(t)}{({v}_{ji}+{v}_{ij}+ { a}_{ij} {M}_{ij}\left({B}_{H}( pH(t))\right){B}_{j}(t))}$$
    (6a)
    $${Q}_{ij}\left({B}_{i}(t),{B}_{j}(t)\right)= \frac{{a}_{ij}{v}_{ij}{{M}_{ij}\left({B}_{H}( pH(t))\right) B}_{i}(t){B}_{j}(t)}{({v}_{ji}+{v}_{ij}{M}_{ij} \left({B}_{H}( pH(t))\right)+{a}_{ij}{B}_{j}(t))}$$
    (6b)

In the EwE package, case (a) is parameterized by making use of the so-called forcing functions (\({S}_{ij})\), which are used to introduce the effects of external variables on trophic interactions and production. Conversely, case (b) is parameterized by making use of the so-called mediation functions (\({M}_{ij})\), which typically simulate the influence of a third (mediating) variable (in this case, the HFS biomass) on predator–prey interactions among the other two functional groups (Christensen and Walters 2004).

Figure 1
figure 1

Diagrams of the food web. Schematic diagrams with the idealized food web model on the left (in black). A change in HFS growth and state (gray arrow) has a direct effect on trophic relationships in which the HFS are directly involved as prey or predator (continuous blue line in A and B) but can also have a non-trophic (mediated) effect on other trophic relationships, not directly involving HFS, for example, refugees for prey and aggregation points for predation ('continuous waved orange line, B). In both cases in due time, the alterations propagate through the food web (dashed line). B Illustrates the cumulated trophic and non trophic effects.

Full size image

The mediation functions can be used to represent a non-trophic effect between species within a food web modeling framework (Figure 1b). For example, mediation functions have been applied to a system in which changes in eelgrass coverage affect the vulnerability of prey to predators through refuge provisioning and/or improving search efficiency by predators through the aggregation of prey (Plummer and others 2013). Moreover, mediation functions have been used to represent the decrease in the vulnerability of juvenile crabs to predators in habitats with increased vegetation (Ma and others 2010), the increases in the feeding area and food availability related to changes in kelp forests (Espinosa-Romero and others 2011), the loss of refuges and consequent negative impacts for juvenile fishes related to the removal of complex substrates.

In detail, three mediation functions were applied to specific functional groups (Table S2) to simulate the loss of habitat services due to the decrease in habitat-forming functional groups (seagrass and corals). These functions represent three types of non-trophic effects within the food web model due to HFS decreases. The first function reduces the number of refuges for small prey species and increases their vulnerability to predators (Figure S3; M1). The second function decreases the food search rate of juveniles due to the loss of prey aggregation (Figure S3; M2), and the third function decreases the search efficiency of predator species linked to the two habitats analyzed (that is, seagrass and coralligenous habitat; Figure S3; M1). The initial mediation value was set to 1 for M1 and M2 and 0.5 for M3 (Figure S3; M3). According to Harvey (2014), the last value was chosen to dampen the strength of the mediation effects on functional groups that are not strictly associated with the HFSs. For example, the effects of HFS reduction might be lower for pelagic predators that can change their feeding grounds by moving to other sites.

Scenarios of OA

We designed a set of 5 scenarios (FF1–FF5) representing the effects of 5 levels of OA on altering the habitat-forming species (HFSs) biomass based on properly defined ‘forcing functions’ Sij(pH(t)). Given the still high uncertainty regarding the shape of the functions that describe these responses, a linear relation between Sij and pH(t) was adopted, as highlighted by Busch and others (2013). The slope (sensitivity to pH) was adjusted to produce HFS biomass decreases of 10%, 30%, 50%, 70% and 90% by the end of 2100 under the 5 different scenarios. The two HFSs (seagrass and coralligenous) present clearly different slopes that have been empirically identified through a number of trial-and-error simulation runs (see Table 1 for a summary of abbreviations).

Table 1 List of Abbreviations.
Full size table

The first set of 5 scenarios described in case (a) did not include any mediation effect and helped us to analyze the dynamic modification of trophic interactions involving HFSs and their cascading effects on the whole food web.

In the second set of 5 scenarios (FM1–FM5), we introduced a mediating effect due to habitat loss (that is, refugee, aggregation of prey) in addition to the ‘forcing function’ defined in the first set of scenarios. To this end, we defined the mediation functions corresponding to the non-trophic effects of seagrass and coralligenous species on other species, as identified in the conceptual model of Zunino and others (2019), based on the specific review of direct and indirect OA impacts on Mediterranean benthic ecosystems (Zunino and others 2017).

Simulations were then run for the same HFS reduction considered previously (scenarios FM1–5, imposed HFS biomass reductions of 10%, 30%, 50%, 70% and 90%, respectively). The comparison of the two sets of simulations allowed us to quantify how the ecosystem and its components responded to changes in the biomass of HFSs and to their direct and mediated effects.

Analysis of the Results of Ecological Network Analysis

Each simulation was run for 100 years, and the results were used to assess the changes in the ecological structure and the degree of maturity of the ecosystem and its functioning under the different scenarios. The selected indicators were total biomass (TB), total system throughput (TST), capacity (C), ascendency (A), relative ascendency (A/C), the proportion of flows to detritus (PFD), the mean trophic level of the community (mTLco), the mean path length (MPL), average mutual information (AMI), ecosystem entropy (H), redundancy (R), the Finn cycling index (FCI), the evenness index (E), the Kempton (Q90) index and total catches (TC) (see Table 2 for indicator descriptions and references).

Table 2 Name, Acronym and Definitions of Indicators Used in the Analysis.
Full size table

The uncertainty in biomasses, catches and indicators were assessed using Monte Carlo (MC) techniques (that is, as the spread of an ensemble of 100 simulations with differing values of most uncertain input parameters, considering a 0.05 coefficient of variation (CV) around the input parameters for biomass (B), production (P/B), consumption rates (Q/B) and ecotrophic efficiency (EE)). The analysis was performed using the EwE Ecosampler plugin (Steenbeek and others 2018), and significant differences between scenarios were assessed using the pairwise Wilcox test in R (R Core Team 2016). Finally, some indicators previously selected in the ENA were used as proxies to provide insight into ecosystem service provisioning.

We estimated the economic value of commercial fishing using market data from the ISMEA database (www.ismea.it), which reports the price of landed species. The market value of each functional group was obtained as the average of the price of each species belonging to the functional group (see Table S1), and we assumed that a change in biomass did not affect harvest costs.

Results

Changes in Ecosystem Biomass

The simultaneous reductions in seagrass and coralligenous HFSs modeled through 100 realizations for each of the 5 scenarios obtained by decreasing their biomass by 10%, 30%, 50%, 70% and 90% resulted in cascading effects that significantly altered ecosystem biomass. These effects were evident on average biomass change for each ecosystem component in both sets of scenarios (FF, without mediating effects, and FM, with mediating effects) (Figure 2, Table S3).

Figure 2
figure 2

Relative changes in ecosystem functional group biomass. Boxplots showing the relative changes in biomass (median and 25th and 75th percentiles) were generated by applying forcing functions (FF scenarios-blue) or mediating functions (FM scenarios-orange) for all scenarios (x-axes).

Full size image

For instance, a reduction of 30% in the habitat-forming functional groups resulted in overall reductions of 28% and 23% in the total biomass of the food web in scenarios FF2 and FM2, respectively, with significant decreases (p-value < 0.001) in many different groups. For example, in FF2, the reduction in HFSs had negative trophic effects on S. salpa as both juveniles (− 18.5%) and adults (− 31%), decapods+ (− 18%), macrofauna+ (− 19%), seaworms (− 10%) and sea urchins (− 22%), which fed directly and primarily on the seagrass meadows. The decreased biomass of the above species led to negative effects on other functional groups for which prey biomass was decreased. Hence, dusky groupers and wrasses+ decreased by 22% and 20% each, scorpionfishes (large scorpionfishes+ and scorpionfishes & combers+) by 15%, pagellus by 16%, striped red mullets+ by 14%, and echinoderms+ by 13% (Figure 2, Table S3). The biomass reduction for the above-mentioned functional groups was even higher in the scenarios with a higher reduction rate of HFSs (FF3, FF4, FF5). The direction of the changes in functional group biomass was consistent across the 5 scenarios for the majority of the functional groups. The magnitude of the change increased linearly with the reduction in HFSs except for a slight increase in the biomass of some pelagic organisms that increased in all the scenarios due to increased zooplankton biomass. The increase in zooplankton biomass was a trophic consequence of the reduction of the main zooplankton predators, the polyps of the gorgonians. For example, in FF2, this increment led to an increase in small tuna (10%), horse mackerel (7%) and sand smelt (5%) (Figure 2, Table S3).

Under the FM scenarios, direct and mediated effects affected many functional groups in the same direction, inducing cumulative impacts that were greater than those in the corresponding FF scenarios. For example, under the FM2 scenario, a 30% HFS biomass reduction induced a greater biomass decrease (− 23% of the total biomass) than under the FF2 scenario (− 28% of the total biomass) in dusky groupers (mean of multistanza groups − 41%), striped red mullets+ (− 18%), pagellus (− 30%), salema (mean of multistanza groups − 37%), decapods+ (− 31%), cephalopods (− 11%), macrofauna+ (− 29%), sea urchins (− 30%), meiofauna (− 22%) and detritus (− 5%). In other cases, mediation and direct impacts acted in opposite directions due to the onset of top-down and bottom-up interactions. For example, in the FM2 scenario, the trophic and non-trophic effects led to positive effects on some functional groups, such as zooplankton (zooplankton-large +25%; zooplankton-small +6%), by increasing their prey biomass, and they positively affected other functional groups by reducing predation pressure (that is, the biomass of gobies+, mullets+ and suspensivores+ increased) (Figure 2, Table S3).

As observed under the FF scenarios, the changes in biomass observed for the functional groups were even higher in the scenarios with a higher reduction rate of HFSs (FM3, FM4, FM5), showing 38% (p-value < 0.001), 53% (p-value < 0.001) and 64% reductions of the overall biomass (p-value < 0.001) (Supplementary material Table S3).

Notably, in both the FF and FM simulations, the variance associated with the change in biomass increased as the reduction of HFSs increased (see Figure 2 and Table S3).

Changes in Ecological Indicators of Ecosystem Status

The scenario analysis detected significant changes in most of the ecological network analysis (ENA) indicators and ecosystem descriptors (p-value < 0.05, see Table S4, Figure 3). No significant differences were found between the scenarios with a 10% reduction of HFS (FF1 vs FM1), except for the evenness index, TLC and MTLco (p-value > 0.05), whereas a slightly significant difference was found when comparing the two scenarios with the original food web—Scenario 0 (p-value < 0.05).

Figure 3
figure 3

Ecological indicators analysis. Boxplots showing the relative changes in indicators values (median, 25th and 75th percentiles). The estimated indicators for reference (gray box plot) and for the scenarios involving HFS reduction by 10%, 30%, 50%, 70% and 90% due to OA. Scenarios were run either assuming only ecosystem trophic effects (FF scenarios, blue) or including HFS non-trophic effects (FM scenarios, orange).

Full size image

The increase in OA implied a decline in the total system throughput (TST), suggestive of a decline in ecosystem size, in all the scenarios showing a significant difference in comparison with the unperturbed condition for all cases except for the two 10% scenarios (p-value > 0.001). No significant differences were found between the FF and FM scenarios for the − 10%, − 30% and − 50% HFS decreases (p-value > 0.001; Table S4, Figure 3).

The development capacity (C) and ascendancy (A) indicators presented a similar trend, which confirmed decrements in ecosystem development and growth potential as perturbation increased. However, a decreasing relative ascendancy (A/C) suggested a slight reduction in the organization of the food web in the FF scenarios only for an HFS reduction up to − 70% (FF4; W = 8314, p-value < 0.001). These effects changed in the FM4 and FM5 scenarios, showing recovery with higher perturbation, possibly due to the increased relative biomass of species belonging to high trophic levels (Table S4, Figure 3).

Indeed, compared to unperturbed conditions, the simulations indicated a significant decrease in the flows to detritus (PFD) in all scenarios, which was reflected by the increase in the mean trophic level of the community (mTLco), indicating a shift of the food web to a higher trophic level (> 2) and a decrease in the relative biomass through detritivorous and herbivorous functional groups. For large perturbations, as in the FF5 scenario, a considerably different value for the mTLco index indicated an increase in the relative biomass of high-trophic-level predators.

Accordingly, the increase in the mean path length (MPL) indicated a change in the system flow structure, with a greater relative abundance at high trophic levels than in the unperturbed system. The results were also consistent with the significant reduction in the average mutual information (AMI) indicator, which indicated that a larger number of groups (high trophic levels) were involved in the transfer of energy and that flows were less constrained to channel energy through specific pathways. This trend was confirmed by the decrease in entropy (H) and the higher redundancy (R) shown in all scenarios compared to the reference, which indicated that the system presented high “strength in reserve” and buffered the effects of perturbation, at least until the − 70% reduction in habitat-forming species.

The results of the structure indicators (TST, A, C, A/C, H) were reflected in the Finn cycling index (FCI), which started to increase significantly under both FF-30% (W = 2544, p-value < 0.001) and FM-30% (W = 2574, p-value < 0.005) (mean +6%) and continued to increase significantly under all the other scenarios. Increased cycling indicates a longer residence time of energy in the system and the onset of some self-reliance of the internal flow organization (Ulanowicz 1983) to cope with the perturbation. The only exception was observed for FCI under extreme perturbation (FM-90%), where FCI decreased by 17% (W = 9090, p-value < 0.001). In fact, under such extreme perturbation, the change in the biomass distribution led to a collapse of the trophic web and its functioning.

Changes in Biodiversity Indicators

The biodiversity indicators of the evenness index, E, and the Kempton Q90 index showed opposite results. Indeed, evenness significantly increased in all scenarios (p-value < 0.001), whereas Q90, a measure of evenness and richness, decreased (except under scenario 1, − 10%; W = 5693, p-value  = 0.09). Such disagreement between the metrics implied that evenness was increased due to the new distribution of energy among many functional groups (and, thus, increased E), while the number of functional groups with biomass within the 5th and 95th percentiles (that is, the proxy for system richness as measured by Q90) decreased more dramatically.

Changes in Catches and their Economic Value

Under the assumption that SSF and the recreational fishing effort were maintained at a constant level, the total catches (TC) significantly decreased in all the perturbed scenarios of OA impacts (p-value < 0.001). The lowest average values were − 29% TC, under FF5 (reduction of HFS by − 90%; W = 10,000, p-value < 0.001), and − 20% TC, under FM4 (W = 9829, p-value < 0.001). In FM5 (HFS reduction by − 90%), the impact was significantly negative (W = 8133, p-value < 0.001), but the associated deviation made the result less reliable. In FF4, for the 16 functional groups listed by Prato and others (2016), all groups except for small tunas+, horse mackerels+ and sand smelts+ responded negatively to the decrease in P. oceanica and coral biomass. Instead, the inclusion of non-trophic effects (FM scenarios) resulted in critical changes in the biomass of some functional groups, sometimes dampening the impacts of the proposed scenarios of OA. In particular, in FM4 (trophic and non-trophic effects of a − 70% HFS biomass reduction), all groups except for amberjack and dentex+, scorpionfishes & combers+, horse mackerels+, sand smelts+, diplodus+, gobies+ and mullets responded negatively. The main reason was the compensatory effects of the FM scenarios, in which the reduced biomass of predators led to the release of meso-predators from trophic pressure that was not highlighted in Scenario FF. In addition to the TC indicator, we have assessed the changes in landing in monetary terms by applying a market-based indicator (the average price of each species in the functional group: ISMEA database) to the commercially important species caught by SSF. We measured a total annual landing value of commercial species within the MPA of €366.5 ha/year (Prato and others 2016) for the present-state scenario based on current market value (Table 3). The SSF is subjected to a loss in economic revenue that is similar in the two scenarios. For example, FF4 and FM4 are both subjected to a reduction of the SSF economic income by 23% (V = 58, p-value = 0.76).

Table 3 Monetary Value Associated with the Commercially Important Species.
Full size table

The market-based indicator (the average price of each species belonging to the functional group: ISMEA database) suggested that the degradation cost associated with OA impacts, as represented in scenario FF4 and scenario FM4, was almost − 20% of the economic income of SSF (Table 3, Figure S4). SSF is of great importance in terms of job opportunities and the contribution to the economies of coastal communities in Europe, generating approximately half of all direct employment within the EU fishing sector and representing a quarter of the catch value (FAO 2018).

Discussion

Our study highlights how OA-induced changes in HFSs such as P. oceanica and coralligenous habitat could have substantial ecological consequences due to complex food web interactions and how these changes adversely impact ecological indicators and fisheries catches chosen as proxies for related ecosystem services. We explored the potential effects of habitat loss due to OA by implementing a new parametrization within a mathematical food web model and describing the decrease in HFS biomass, the cascading potential trophic and non-trophic effects on other species, the related consequences in terms of ecosystem structure and functioning, and eventually the loss of ecosystem services. The modeling approaches can help the understanding and quantification of ecosystem trends in future scenarios, providing indications that could help both prevention and mitigation of climate change impacts.

Ecosystem Effects of OA through HFS

Our results show that the habitat changes due to OA could threaten keystone species (that is, dusky groupers, amberjack & dentex+ and large-scale scorpionfish+, Prato and others 2016) by reducing both their foraging area and their habitat, thus decreasing the services that these habitats provide to these species (that is, the protection offered in nursery areas). A slight decrease in the biomass of medium and large dusky groupers will lead to large indirect effects on the food web that should favor cephalopods+, which are the groupers’ favorite prey.

Furthermore, a decrease in gorgonians will lead to an increase in the biomass of their zooplankton prey. In turn, the increase in zooplankton will contribute to an increase in the biomass of sand smelts+, showing an unexpected strength of the cascading effects between these trophic links, which are much stronger than the non-trophic effects imposed. It is worth noting that these results are often counterintuitive due to the complex interactions between trophic relationships impacted in a bottom-up or top-down manner. It was not possible to consider the responses of all functional groups to OA. However, major trophic effects induced by the reduction of HFS biomass were considered, and our results highlight the importance of addressing the anthropogenic impacts in complex systems by considering the whole food web. Indeed, competition and predation relationships are able to alter the community responses to habitat degradation due to OA in complex ways. Nevertheless, we must note that this model dynamically accounts for the potential adaptive shifts in the diet (see also Walters and Christensen 2007) of the impacted functional groups only partially and that a high percentage of production is exported outside the system. For instance, detritus accumulation due to HFS degradation might be underestimated, and all detritus-based flows may be insensitive to functional group changes.

Inclusion of Habitat-Mediated Effects

Contrary to what we expected, the scenarios representing both trophic effects (FF) and trophic and non-trophic effects (FM) indicated a decrease in cephalopods induced by a decrease in their prey (that is, decapods+ and macrofauna+) that are directly dependent on P. oceanica. The reduction in HFSs (FF) simulated together with the introduction of non-trophic effects (FM) indicated that critical changes in the biomass of some functional groups sometimes dampen the impacts of the proposed scenarios of OA. This occurred because FM scenarios affecting predator–prey relationships led to a greater reduction in HFS predators (that is, sea urchins, salema, decapods+ and especially diplodus+), which in turn favored some recovery of the HFSs. The analysis of biomass changes under the FF and FM scenarios highlighted the synergistic effects of habitat-mediated and trophic OA effects on horse mackerel, decapods, dusky grouper, meiofauna, pagellus and salema. These groups showed increased effects when habitat mediation was included. Conversely, habitat-mediated effects acted antagonistically to OA effects for sea worms, scorpionfish and scombers, wrasses and striped red mullet; for these groups of species, the habitat effects were beneficial for reducing OA impacts. Complex habitats enhance the aggregation of prey, and their loss may lead to a decrease in the search rate of predators (that is, loss of habitat aggregation points for prey). On the other hand, these complex structured habitats offer protection to fishes and invertebrates in both the juvenile and adult stages, and the loss of these functions may increase the vulnerability of some species to predation (that is, loss of prey refuges). It is documented that predation upon juvenile or invertebrate benthic species is significantly lower in seagrass than in unvegetated habitats (Espinosa-Romero and others 2011). However, since our model represents only a few ontogenetic phases, we are not able to identify all of the different effects of trophic and non-trophic impacts caused by the loss of HFSs on the functional groups in the model. The underrepresentation of non-trophic effects is counterbalanced by the simplistic assumptions made for them, which may lead to an overestimation of the effects of habitat loss.

Effects of OA on Ecosystem Development Indicators

The above-reported changes in the examined indicators suggest that under the effects of pH increases and reductions in HFS biomass, the system decreases in size (declines of TST, C and A) but also reorganizes itself in such a way that relatively higher biomasses occur at higher TLs (increase in mTLco), to maximize the use of energy by increasing energy passage (increased MPL) and increasing ecosystem redundancy and strength overall (increase of R). As a result, energy cycling in the food web is increasing (FCI increasing). Thus, overall, the system reorganizes itself under perturbation.

Despite the growing empirical evidence of OA effects on marine ecosystems and their component species, the prediction of the impacts of OA across the entire ecosystem remains challenging. Indicators that integrate all the species-specific responses at local scales and the cascading effects on the entire ecosystem support the view that highly structured ecosystems such as that of the Portofino area exhibit high levels of resilience to OA perturbations. Indeed, the results show a system reorganization under OA impacts that partially counterbalances individual species degradation and modification. Species habitat affinities, physiological preferences, life histories and interspecific interactions vary among systems and are rarely known in sufficient detail to make predictions, but the overall signal at the ecosystem level in terms of development indicators clearly reveals the great potential for such systems to respond and readjust under OA. Additional analyses are needed to evaluate whether such ecosystem compensatory effects might be less pronounced in areas that are not protected. However, counterbalancing the perturbations by maintaining a high system development status sensu Odum (1969) does not necessarily mean that all processes of interest to humankind are preserved.

Ecosystem Services Under OA

The analysis of the set of proxies used to assess changes in ecosystem capabilities to provide ESs showed that future OA-impacted ecosystems will present a lower capacity to provide provisioning, cultural, support and regulatory services.

In particular, we used the Kempton Q90 biodiversity indicator as a proxy for impacts on supporting services (according to the MEA 2005 classification), resulting in mild negative impact under both (FF and FM) the − 30% and − 50% scenarios.

We selected structural indicators (that is, TST and ascendency) as proxies of regulatory services since they provide information related to the extent of nutrient cycling and biological control over these services (Hattam and others 2015). These services also appear to be reduced under increased perturbation.

The fraction of the total biomass related to commercial fishing (the potential fished biomass) is used as a proxy for provisioning services, and it is moderately impacted under all the perturbed scenarios (Figure 4).

Figure 4
figure 4

Ecosystem services risk impact: the qualitative assessment of the OA impacts on the provisioning of the services (y-axis) derived from the quantitative responses of the selected indicators. The strength of the impact was categorized according to the legend. The services that are expected to be highly impacted/lost under the degradation of the HFS are indicated in orange (that is, supporting services). Yellow indicates the services that are expected to be moderately reduced (provisioning, regulating and cultural), whereas the undetectable changes are indicated in white.

Full size image

Finally, the selected indicators of cultural services were the fraction of total catches, which is related to small-scale and recreational fisheries, since they have high cultural, historical and socioeconomic significance (Prato and others 2016), and the biomass of HFS, which was identified as an important cultural parameter in a study on divers preferences (Zunino and others 2020) (Figure 4). Model projections indicate decreases in these proxies.

Figure 4 highlights that all ESs decrease with increasing acidification.

Conclusion

Climate change impacts together with other anthropogenic pressures (that is, fishing, pollution, and eutrophication) are threatening marine ecosystems, decreasing their resilience and their ability to provide benefits and well-being to human society. Despite the widely recognized functions and services provided by coralligenous species and P. oceanica, few studies have attempted to quantify the value of these habitats in supporting ecosystem services.

The analysis of scenarios of increasing OA through 2100 and its impacts on HFS clearly highlighted increasing impacts on food web components at higher OA levels, in which the effects of HFS reduction propagate through the food web with significant negative cascading impacts. Trophic and non-trophic impacts of HFS reductions showed synergistic effects, suggesting unexpected very high overall effects that might be difficult to predict.

The evaluation of such changes through ecological network indicators reveals that HFS-structured ecosystems undergo a reorganization of their functioning by decreasing in size but maintaining good developmental stages (sensu Odum 1969), high ecosystem organization, high average trophic level of the community and high average energy passages. This suggests that healthy ecological systems are able to cope with and adapt to some changes in OA to maintain their ecological functioning. Conversely, proxies for ecosystem services clearly indicate a significant decline due to OA impacts on HFSs in the long term. The most relevant impacts of OA on the ecosystem are therefore not necessarily related to its health and stability but take the form of a reduced ability to provide ESs. These results converge to highlight the crucial role of HFSs in marine ecosystem functioning, especially in demonstrating that some OA effects might cause alterations in marine ecosystems whose effects are negative mainly from the perspective of humankind and therefore supporting the conclusion that the reaction of “nature” to OA needs to be accompanied by a further adaptation by “humans.”

Although these findings were obtained in an area in the northwestern Mediterranean Sea, they are likely to show more general validity.