Impact of Ocean Acidification on Ecosystem Functioning and Services in Habitat-Forming Species and Marine Ecosystems

Ocean acidification (OA) is expected to impact habitat-forming species (HFS), with cascading effects on the whole marine ecosystem and related services that are seldom quantified. Here, the changes in HFSs biomass due to OA are modeled using a food web ecosystem model, and the trophic and non-trophic cascading effects on the marine community are investigated. The food web model represents a well-studied coastal marine protected area in the NW Mediterranean Sea where coralligenous reefs and Posidonia oceanica meadows constitute important HFS. The model is used to implement 5 scenarios of habitat degradation, that is, reduction of HFS biomass, induced by increasing OA and to quantify the potential changes in ecosystem properties and indicators of ecosystem services over the next 100 years. The changes in ecosystem indicators highlight a decrease in the size of the system and a reorganization of energy flows suggesting a high degree of ecosystem development. All the proxies for ecosystem services show significant decreases in their values. Although representing only a portion of the possible impacts of OA, the findings are consistent with the idea that ecological systems can react to OA effects to maintain the level of ecosystem development, but the new organization might not be optimal from an anthropocentric viewpoint.

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 nontrophic effects (hereafter called ''mediating effects'') on the biomass of ecosystem components, ecosystem processes, interactions, and properties. Moreover, this work quantifies the potential longterm 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 wholeecosystem 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 smallscale 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 (B i ) over time (t) as follows: where dB i /dt represents the change in the biomass of group i occurring over a time interval of dt; g is growth efficiency; Q ji defines the trophic relationships between a consumer (j) and a prey (i); I i and e i are the immigration rate and emigration rate, respectively; F i is the fishing mortality rate; and M i summarizes all other mortality not explicitly represented in the model by predation and fishing mortality. PP i is the function that drives the primary production of autotrophs over time, as follows: where r i is the maximum production rate, and k i 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 B i (t), the predator biomass B j (t), the effective search rate for species i by species j, a ij , and the rates at which species i moves between the vulnerable and invulnerable prey pools, v ij (Christensen and Walters 2004): 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: (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 r i in Eq.
(2) with a function, r i S(pH(t)), which depends on the timevarying value of pH(t): Similarly, for heterotrophs, the constant parameter a ij in Eq. (3) is replaced with a function that depends on pH(t), a i,j *S ij (pH(t)): OA Impact on Ecosystem Functioning and Services (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 a ij , describing the selected trophic relationships between species i and j (different from h), in Eq.
(3) with the function a ij AEM ij (B h (pH(t)), which depends on the biomass of the habitatforming species h. Analogously, the constant parameter v ij is replaced by the function v ij AE-M ij (B h (pH(t)): 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). 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 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. 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' S ij (pH(t)). Given the still high uncertainty regarding the shape of the functions that describe these responses, a linear relation between S ij 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).
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).
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 The sum of all flows in a food web, representing a measure of the ''size'' of a system in terms of flow''. T ij is the flow between two compartments Finn (1976) Development capacity A measure of the potential of an ecosystem to develop and the theoretical maximum ascendency (A) Describes the growth and development of the system. T i is the sum of all material leaving the ith compartment; T j is the sum of all flows entering the jth compartment Ulanowicz (1986) Relative ascendency A/C Index of organization of the food web Ulanowicz (1986) Redundancy or internal flow overhead Indicates the change in the degrees of freedom of the system or the distribution of energy flow pathways in the system. It is considered a measure of ecosystem pathways

TST
Represents the diversity of flows in a system, providing an indication of the total uncertainty embodied in the given configuration of flows Ulanowicz (2004) Mean path length MPL = TST P inputs Accounts for the number of functional groups involved in a flow of matter and represents the average number of groups that an inflow or outflow passes through Finn (1976) Kempton Q-90 index Measures species diversity in the whole ecosystem. This statistic represents the slope of the cumulative species abundance curve between the 10th and 90th percentiles. Each functional group represents one species, and the biomass is a proxy for the number of individuals of that species. (S) is the total number of functional groups in the model; R 1 and R 2 are the representative biomass values of the 10th and 90th percentiles in the cumulative abundance distribution) Ainsworth and Pitcher (2006) Evenness E ¼ À P n i p i lnðp i Þ Shannon's index is widely used in ecology as a measure of species diversity. Since the number of functional groups is fixed, the index applied to food web results essentially measures evenness Ulanowicz (1986) Finn's cycling index FCI ¼ TSTc TSTtot Calculates the proportion of the total system throughput (TST) that is recycled into the system Finn (1976) Mean trophic level of the community MTLco Weighted average trophic level for functional groups with TL > 2, which is used to reduce the variability in terms of biomass and production that comes with using low trophic levels Pinnegar and others (2002) Total catches RYj Sum of all fishery catches extracted from the modeled ecosystem in a given year 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.

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).
For instance, a reduction of 30% in the habitatforming 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 & com-bers+) by 15%, pagellus by 16%, striped red mul-lets+ 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).
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).
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) 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). 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-trophiclevel 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 (pvalue < 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 nontrophic 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 mack-erels+, 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 marketbased 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 com-mercial species within the MPA of e366.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).
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 nontrophic 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 diplo-dus+), 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. Indica-tors 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).
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. FUNDING Open access funding provided by the Acid.it Project ''Science for mitigation and adaptation policies to ecological and socioeconomic impacts of acidification on Italian seas''.

OP E N A C C ES S
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copy-right holder. To view a copy of this licence, visit h ttp://creativecommons.org/licenses/by/4.0/.

D A T A AV AI L ABI LI TY
The original EwE model developed by Prato and others (2016) is uploaded on the Ecobase database. The results of this study are available in the supplementary material. 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.