Exploration of multiple post-extinction compensatory scenarios improves the likelihood of determining the most realistic ecosystem future

A research agenda is currently developing around predicting the functional response of ecosystems to local alterations of biodiversity associated with anthropogenic activity, but existing conceptual and empirical frameworks do not serve this area well as most lack ecological realism. Here, in order to advance credible projections of future ecosystems, we use a trait-based model for marine benthic communities to inform how increasing trawling pressure changes the biological-mediation of seabed functioning. Our simulations reveal that local loss of species, and the associated compensatory community response, lead to multiple and disparate biogeochemical alterations that are contingent on relative vulnerabilities to extinction, environmental and biological context, and the level of functional redundancy within replacement taxa. Consequently, we find that small changes in faunal mediation caused by community re-organisation can disproportionately affect some biogeochemical components (macronutrients), whilst having less effect on others (carbon, pigments). Our observations indicate that the vulnerability of communities to future human-induced change is better established by identifying the relative magnitude and direction of covariance between community response and effect traits. Hence, projections that primarily focus on the most common or most productive species are unlikely to prove reliable in identifying the most likely ecological outcome necessary to support management strategies.


Introduction
The widespread reorganisation of ecological communities associated with anthropogenic activity and climate change are heightening concerns about the likely consequences for ecosystems [1,2]. Despite some consistency in the evidence over the directional effects of biodiversity on ecosystem functioning [3], a disconnect between field observations and expectation based on theory and empirical investigation [4] frustrates efforts to project the most likely ecosystem futures. One major difficulty is that community responses to disturbance are not random and reflect a variety of compensatory interactions amongst surviving species that are dependent on prevailing conditions [5]. Hence, the response of a community is largely mediated by individual response traits that define a species capacity to withstand net environmental conditions [6] and, to a lesser extent, by their effect traits that govern how organisms interact with the environment [7]. Recent simulations suggest that the most likely ecological outcome is conditional on the collinearity amongst response and effect traits [8] and how these traits covary with the risk of extinction imposed by short-and long-term environmental forcing [9]. This means that inherent differences in environmental history [10,11], as well as species-specific variation in functional contribution [12], level of extinction risk [13] and capacity to compensate [14][15][16] decouple the response of the surviving community from the level of taxonomic reorganisation [17,18]. Hence, it is important to recognise that the capacity of species and communities to respond to environmental influence can be constrained by the trophic structure and level of functional redundancy, in particular the degree of positive versus negative species interactions [8,19], that are expressed within the surviving community. These, in turn, are dependent on environmental context [20,21], biological circumstance (e.g. competitive release) [22,23], and/or the a priori experience of species to recent abiotic and biotic conditions [24] that can lead to individual changes in behaviour [25]; all of which are compounded, and can change in relative importance, at larger temporal and spatial scales [26,27]. Accounting for these sources of variation when considering community functional responses to extant or anticipated perturbation have, however, received little attention despite calls for the development of more realistic scenarios [28] that are grounded by theory and/or field observation [29,30]. A more systematic integration of these features in conceptual and empirical research frameworks is needed, because it is likely to accelerate the implementation of practical solutions [1].
Recent developments in theory [14,31,32], observations [33] and emerging experimental evidence [15,16] emphasise the importance of population dynamics and compensatory mechanisms in moderating the ecological repercussions of altered biodiversity, but these mechanisms have seldom been evaluated in natural systems undergoing progressive forcing [34]. Intuitively, the vulnerability of a natural system is a function of the type and the extent of the response of the surviving community which, in turn, is ultimately contingent on the perturbation that has triggered the response and the wider environmental context [24,35]. Considering appropriate compensatory features within a projection, therefore, is difficult to achieve in practice because the response of the community is not always foreseeable; it is controlled by a complex array of interwoven drivers that are not yet fully understood and which may include some level of ecological surprise [36,37]. Common practice has been to explore a range of possibilities that span all imaginable outcomes [14] but, as a means of homing in on the most likely outcome, this approach lacks an objective way of scenario selection and, in any case, insufficiently captures the complexity and context-dependence of an ecological response. Investigations of species performance along natural diversity gradients [28,29] have been instrumental in narrowing the number of possibilities, and are useful in stimulating reflective debate about what constitutes a most likely scenario [38]. From these narratives, refined and refocussed experimental designs have emerged which are beginning to provide results more aligned with realistic biodiversity loss [39][40][41]. Impetus to formulate better model projections [42,43] that encapsulate, preserve and propagate this realism from small-scale experimentation to the landscape level of natural systems [44,45] is needed to support decision makers and policy platforms [46].
Expert knowledge is a critical factor in improving the ecological realism of modelling frameworks and tools [42], but requires the merger of perspectives from different disciplines [47,48]. Recent extensive integrative research programmes that embrace this philosophy (e.g. Shelf sea biogeochemistry; [49]) , allow mechanistic insight (i.e. ecosystem processes) to be combined with environmental and ecological status (i.e. communities), and confront models with empirically grounded information to generate stronger predictive relationships [50].
Here, we use a trait-based model of marine benthic communities to investigate change in the biologicalmediation of seabed functioning (carbon and nutrient cycling) under an increasing gradient of chronic trawling pressure. We control for confounding variables by using empirical data collected under comparable environmental conditions at the same time of the year, whilst allowing for the influence of context by constraining the model by sediment type. Our a priori expectation was that the biological mediation of carbon and nutrient cycling would differ with sediment type, irrespective of species composition, and that projections of the ecosystem consequences of altered biodiversity within these environments would not follow the same trajectory pattern. Further, we anticipated that the average body-size of the community would decrease with increased levels of trawling pressure [51,52] and, because of the importance of body size in determining a species functional contribution [53,54], would lead to an accelerating decline in seabed functions as trawling pressure increases. Realisation of these expectations would emphasize the importance of context in determining the ecological consequences of altered biodiversity. In addition, it would highlight the need to identify and incorporate sources of post-extinction compensatory dynamics and environmental change in projections of biodiversity-ecosystem futures.

Study location and model parameters
We use a comprehensive benthic survey of 55 stations in the southern part of the Celtic Sea shelf (∼2200 km 2 , 50°16.37-50°29.62N, 6°8.03-7°11.41W; [49]; data S1 and table S1 (available online at stacks.iop.org/ERC/3/ 045001/mmedia)), to parameterize models that predict how alterations to biodiversity associated with trawling pressure affect seabed function. The selected stations contrast in sediment type and macrofaunal community structure, but exhibit minimal variation in water depth, temperature, salinity, and tidal/wave dynamics. To account for differences in the level of forcing a community may experience, we matched each station to estimates of trawling pressure (table S2) derived from aggregated vessel monitoring system (VMS) data for bottom trawled gears (vessel speed 1-6 knots from 2009 to 2014, standardised by year; [55]).

Modelling tool
Using a probabilistic trait-based model developed for exploring the effects of local extinction scenarios and the associated compensatory response of natural communities [9,14,56], we predict how altered diversity associated with trawling activity is likely to affect seabed functioning. We establish the relationships between an index of community-level bioturbation potential (BPc; [9]), estimated from per capita contributions of sediment-dwelling invertebrates to sediment reworking based on average body-size, abundance, activity level (Mi) and sediment reworking mode (Ri), and measurements of representative biogeochemical metrics (24 biogeochemical parameters; tables 1, S1). Statistically significant relationships (mud, n=6; sand, n=8; table 1) were used to explore how changes in trawling pressure most likely alter faunal mediation of biogeochemical processes at the seafloor.
As bottom trawling activity disassembles communities through the selective removal of specific taxa [57] and influences the compensatory response of the surviving community [58,59], our model selectively eliminates taxa from the immediate species pool before simulating the response of the surviving community through compensatory mechanisms established for the regional species pool (achieved using a probability-based order derived from a gradient of trawling pressure; code S1). Specifically, as taxa are sequentially extirpated and the surviving taxa compensate (increase in biomass) to accommodate the release of resources, a revised community bioturbation potential is calculated, and empirically derived relationships are used to estimate associated biogeochemical metrics. The simulation continues until all taxa within the assemblage become locally extinct. As alteration of natural communities associated with trawling pressure may be offset by taxa from a wider area [60], we allow for taxa present in the regional species pool that were not present in the starting assemblage to be introduced and compensate. This added feature explicitly recognises that a natural community can be qualitatively different along a gradient of disturbance as the most vulnerable species are progressively replaced by less vulnerable taxa from the surrounding area [57,61]. Inclusion of taxa that are present (abundance greater than zero and at risk of local extinction) or absent (abundance equal to zero, no risk of local extinction) allows for the possibility that taxa from the absent pool can arrive and contribute to the present pool as niche space is Table 1. The relationship between bioturbation potential of the community (BPc) and biogeochemical metrics in mud and sand sediments based on a spatial survey of 55 stations in the Celtic Sea (see tables S1-S2). The direction of change with increasing BPc (slope), type of transformation, variance explained (R 2 ) and statistical significance (p-value: * >0.05, ** >0.01 and *** >0.001) are indicated.  released. Once a taxon becomes locally extinct, however, our models do not allow the extirpated taxa to compensate a second time as we assume conditions are no longer supportive (code S2).

Simulation scenarios
As the impact of trawling pressure will depend on local community structure, and because biogeochemical processes are strongly influenced by local environmental conditions [62], we categorised stations by sediment type (mud versus sand). We estimated the consequences of trawling pressure on sediment-dwelling invertebrate communities by deriving a vulnerability for each taxon based on proportional representation across the gradient of trawling pressure. For each sediment type, we defined four levels of trawling intensity (F1, 100-300 h.kw/y, hereafter low; F2, 300-1300 h.kw/y, hereafter medium; F3, 1300-2000 h.kw/y, hereafter moderate; F4, >2000 h.kw/h, hereafter high) and a reference (F0, 0-100 h.kw/y) location (table S2). To determine a taxonspecific vulnerability for the focal level of trawling pressure, we calculated the difference between the mean macrofaunal biomass and abundance of each trawling intensity level within each sediment type and the reference location (higher differences indicate increased levels of local extinction). We used these taxon-specific vulnerabilities to establish the probability of local extinction and, reciprocally, compensation of each taxon. This means that a taxon with a high vulnerability score has both a high probability of going extinct and a low probability to compensate. This methodology was used for each of the 16 trawling scenarios (comprising biomass-and abundance-based vulnerability, two types of sediment and four levels of bottom trawling intensity ;  table S3, table S4), each run for 1000 simulations scenario −1 (code S1).

Statistical analyses
We used a series of linear models to derive the metric-specific relationship between community-level bioturbation potential (BPc) and indicators of biogeochemical performance (n=24) in sand and mud sediments (table S1). Similarly, the vulnerability of taxon-specific contributions to biogeochemical performance was evaluated using linear models to explore the relationships between BPi (bioturbation potential of an individual) and taxon vulnerability (abundanceor biomass-based; table S3). Initial linear models were assessed visually for normality (Q-Q plot), homogeneity of variance (plotted residual versus fitted values) as well as for influential data points (Cook's distance). Data were log or arc sine transformed to improve fit and satisfy parametric statistical assumptions and, where necessary, a constant was added to allow transformation.
To account for the non-linear nature of the response of the biogeochemical projection along the gradient of decreasing species richness, we summarise each simulated biodiversity-function relationship using Generalised Additive Models (GAMs). By using a smoothing function of species richness with trawling pressure as an interaction term, the relationship between biogeochemical metrics and species richness can differ under each intensity of trawling pressure. We also calculated the Coefficient of Variation (CV) around each prediction. This is preferential to other commonly used measurements of dispersion, such as standard deviation, because CV permits the comparison of all of our biogeochemical metrics free from scale effects.
All statistical analysis, data exploration and plotting were performed using the R statistical and programming environment [63] and the R packages 'mgcv' (Generalised Additive Models; [64]) and 'tidyverse' (data exploration and plotting; [65]).

Biological mediation of the seabed biogeochemistry
We find that bioturbation potential relates to a subset of biogeochemistry metrics, dependent on sediment type (mud, 6 significant, 18 non-significant; sand, 8 significant, 16 non-significant; table 1). Only three biogeochemical metrics (phaeopigment, [TOxN]-L, ΔSiO 2 -U) related to bioturbation potential in both mud and sand. In general, we achieved a better goodness of fit in sand (R 2 , range −0.46-0.728, median=0.288) than we did in mud (R 2 , range −0.41-0.426, median=0.027), although it is not possible to make direct comparisons because of differences in the number of stations between sediment type categories (code S1). An increase in BPc was associated with a decrease in particulate carbon-related metrics in mud (Linear regression slope: OC, −0.17; ON, −0.02; Phaeo, −0.43; table 1) but an increase in particulate carbon-related metrics in sand (Linear regression slope: Chla, 0.50; Phaeo, 0.59; table 1). In addition, higher BPc values in both mud and sand showed a tendency to be associated with a greater slope (i.e. greater rate of change in concentration) for silicates in the uppermost region of the sediment profile (ΔSiO 2 -U) and increase in SiO 2 concentration in sand ([SiO 2 ]-U and [SiO 2 ]-L), suggesting that silicates are being relocated from depth to shallower layers. The increase in the slope for ammonium (ΔNH 4 -U) may also indicate a similar transport mechanism. BPc is negatively correlated to total oxidised nitrogen (TOxN) in the lower region of the sediment profile ([TOxN]-L) in both mud and sand sediments and, as the slope of nitrate (ΔNO 3 -U and ΔNO 3 -L) decreases across the sediment profile in sand (although not in mud), it suggests that increased rates of denitrification are occurring closer to the sedimentwater interface. We found no evidence that phosphate-related metrics (PO 4 ) are associated with changes in BPc (table 1).

Vulnerability of taxa
Abundance-based taxa vulnerability indicates, for the most part, that taxa-specific contributions to BPc are not expected to change with trawling pressure (insignificant regression slopes, figure 1). However, when vulnerability was estimated using biomass, we find that taxa that contribute least to bioturbation are selectively removed whilst those that contribute most are more likely to remain (significant negative regression slopes, figure 1). We provide the full list of BPi values ordered by vulnerability to trawling pressure for all 16 local extinction scenarios (figure S1).

Projected ecosystem futures
Our models predict that the loss of taxa and concomitant compensatory community responses associated with changes in trawling pressure lead to substantive changes in multiple biogeochemical metrics in both mud (figure 2, for full probabilistic distribution, see figure S2) and sand (figure 3, for full probabilistic distribution, see figure S3) environments. Irrespective of scenario, most of the change in biogeochemical performance occurs as the first few taxa (typically less than 10) are extirpated, with further species loss having very little additional effect on the level of biogeochemical functioning (solid line, figures 2 and 3). Coefficient of variation of the probabilistic distribution for each level of species richness indicates that these biogeochemical outcomes are not subject to substantial variation other than at the lowest levels of species richness (shaded areas, figures 2 and 3). Hence, based on our model projections, the expectation is that increased trawling pressure will lead to a decrease in mud, and increase in sand, of carbon-related metrics, a decrease in total oxidised nitrogen in mud and sand, a decrease in nitrates in sand, and an increase in ammonium in mud. Silicate concentrations are projected to increase across the whole sediment profile in sand and in the uppermost layers of mud. However, it is notable that some biogeochemical metric projections are less variable than others (lower CV values, figures 2 and 3). In mud, for example, carbon and pigment-related metrics exhibit consistently lower CV values relative to those of the nutrients. Hence, a small change in BPc can translate into a disproportional change in nutrients whilst, at the same time, having little effect on sediment carbon and pigments (table 1, figure S2). In sand, with the exception of the silicates where projections are more variable, most metrics have comparable CV values (table 1, figure S3).
Relative differences in the magnitude of change in trawling pressure are projected to be minimal, both in terms of main trend and predictability. In most scenarios, and for all biogeochemical metrics, outcomes align irrespective of the level of trawling pressure.

Functional contributions of surviving taxa
Our projections across all scenarios indicate that the taxa contributing most to biogeochemical metrics (mud, figures 4; sand, 5) transition from those of the pre-extinction community (in mud, the ribbon worms Nemertea,  the decapod crustacean Goneplax rhomboides, the polychaete Nephtys histricis; in sand, the polychaetes Ditrupa arietina and Glyphohesione klatti), to those of the post-extinction community (in mud, the peracarid crustaceans Ampelisca spinipes, Urothoe elegans and Vaunthompsonia cristata, the polychaetes Scoloplos armiger, Aonides paucibranchiata and the bivalve mollusc Abra nitida; in sand, the polychaetes Spiophanes kroyeri, Prionospio fallax, Scalibregma inflatum and Glycinde nordmanni). Taxa that were important in the transitory phase (in mud, In the upper panels, colour intensity (cold to warm coloration; grey-blue) reflects an increasing density (low to high) of data points, whilst in the lower panels colour shading (low-high, white-dark blue) represents the relative contributions of individual taxa to BPc at each sequential level of local extinction, only a subset of twenty taxa that overall contribute most are represented. In the upper panels, colour intensity (cold to warm coloration; grey-blue) reflects an increasing density (low to high) of data points, whilst in the lower panels colour shading (low-high, white -dark blue) represents the relative contributions of individual taxa to BPc at each sequential level of local extinction, only a subset of twenty taxa that overall contribute most are represented. the polychaetes Harmothoe antilopes, Scalibregma inflatum and the urchin Echinocyamus pusillus; in sand, the polychaetes Dipolydora coeca, Oxydromus flexuosus and the peracarid crustacean Vaunthompsonia cristata) introduce variability in BPc (top panel in figures 4 and 5). The distinction between these transitions is clearer under abundance-based vulnerability than it is under biomass-based vulnerability, where each contributing taxon is better delineated and effects are less variable. The rather flat pattern of BPc values following initial taxa losses indicate high levels of functional redundancy in the intermediate subset of taxa. Nevertheless, a consistent feature of the projected post-extinction communities, regardless of scenario, is that the surviving assemblage yields a higher BPc relative to that of pre-extinction condition, even when replacement taxa have a lower BPi (figure S1).

Discussion
Using a trait-based model for marine benthic communities that supports improved levels of ecological realism, we have demonstrated that species loss associated with changes in trawling pressure can lead to previously unidentified contrasts in the biological-mediation of the seabed that are likely to have significant functional repercussions. Our model projections confirm the role of the surviving community in compensating for the loss of, or reduction in, functionally important species [14,15], and reiterate that the degree of community response reflects the level of covariation between species traits and extinction risk [8,9]. However, our data also emphasise the importance of context [10,20] in determining species vulnerability [34] and post-extinction response trait expression [12]. This means that the reorganisation of communities associated with a specific type or magnitude of disturbance is not pre-determined, but subject to variation, such that the resulting level of functioning depends on differences in the multiple underlying mechanisms that control, for example, the various facets of sediment biogeochemistry. This interpretation is consistent with previous observations in the literature where, despite some aligned outcomes (aRPD; [66,67]; carbon and pigments; [56,68,69]; macronutrient; [70][71][72]; denitrification; [35,69,70,73]), differences in context within [35,68] and across [70,73] spatial areas explain apparent discordance between bioturbation potential and various metrics [69,73], mechanisms [69,74,75] and/or components of seabed function [76]. Lack of clear and consistent relationships between traits may render a predictive approach difficult but such instances are not expected to be common in natural systems [77]. In reality, many components of context moderate the relationship between community reorganisation and ecosystem response, including stressor-related extinction risks [78,79], trophic relationships [80,81], population dynamics [82,83], species diversity [84], resource availability [85] and geography [86]. These contextual actors are often seen as something that must be rigorously controlled for, or excluded, but embracing this source of variation when implementing realistic scenarios of local species loss is likely to be beneficial [39]. Moving forward, the challenge becomes one of identifying which plausible [77,87] and relevant [88] combination of traits are the most appropriate for specific ecosystem functions [69,70], whilst recognising that metric-function relations might not necessarily be generically relatable. Elected traits are likely to be response- [89] and function- [78] specific and, as our results indicate, the assessment of species vulnerability may demand measures based on abundance and/or biomass.
A consistent feature of our model projections, irrespective of scenario and in common with field observations, is that the loss of species leads to an increase in functioning [35,59]. A naïve interpretation might be that this is counterintuitive and contradictory to established consensus [4], but the importance of negative covariance between response and effect revealed in our simulations indicates that the faunal mediation of biogeochemical cycling is not solely governed by positive species interactions [90], it can be further moderated by changing circumstances that alter species composition, dominance and dynamics [12,21,91,92]. Increasing bioturbation activity stimulates (older) organic matter mineralisation in mud sediments, but a transition to interface feeding in sand sediments promotes the translocation of (fresh) carbon from the water column into the sediment profile [93]. Furthermore, shallower macronutrient processes indicate faster metabolization rates that are commonly related to an increased aerobic environment associated with faunal mediation [94] but known to be more critical in mud than in sand where advection prevails [95,96]. Here, in accordance with previous findings [9,53,54], we find that trait biomass strongly associates with functioning such as biogeochemical cycling and that taxa-level biomass is proportional to the relative contribution that each taxon makes to community function [97]. In mud, the compensating community mostly belongs to the same functional group as the extirpated taxa, such that the net gain in function corresponds to a gain in average individual body-mass that follows a sharp decrease in the abundance of a small-sized taxon (ribbon worms, Nemertea, 5 mg) relative to other much larger co-dominant species (e.g. decapod crustacean Goneplax rhomboides, 464.6 mg; gastropod Euspira fusca, 222.6mg) from the starting community. In sand, however, the simulated gain is best explained by transition in functional group dominance, from less mobile species that inhabit the uppermost layers of sediment to species that are more mobile and generate deeper mixing. As the extinction continues beyond the first few extirpations, these effects diminish as the functional contrasts between species are minimised. Hence, the compensatory response of the remaining community yields a higher level of functioning that is largely insensitive to further species loss because of high levels of functional redundancy that buffers against further change in ecosystem properties [34,98].
Mechanisms generating functional redundancy are known to depend on environmental context and the factors driving local extinction which means, in turn, that contrasting organism responses will reflect the driver of species loss [99,100]. Although some generalisation exists, taxa bearing traits related to ecological generalisation and smaller size are expected to fare better under global change [6,101]. Other driver-response relationships can be more specific; for example, temperature rise promotes fast life-histories [102,103], acidification favours the lack of calcified structure [104] and direct development (i.e. non-planktonic) and crawling movements increase as salinity decreases [105]. Our projections suggest that compensation following bottom trawling is primarily undertaken by taxa belonging to the same functional groups (slow biodiffusers, Mi/Ri 3-4; surficial modifiers, Mi/Ri 3-2; and conveyor belt feeders, Mi/Ri, 2-3, sensu [9]). This is broadly in agreement with expectations that sessile or low mobility taxa living at and/or around the sediment-water interface (Mi and Ri lower than 2) are predominantly impacted by bottom trawling [57,106] and is consistent with previous findings [52,57] that emphasise the importance of functional redundancy exercised by smallersized species. Nevertheless, functioning does not appear to be as sensitive to extinction per se as it does to the loss of specific taxonomic groups because high levels of functioning have been reported for low levels of species richness [17,97,107]. Fewer taxa, however, means that the overall functional contribution of the community largely stems from a smaller subset of available trait combinations [104]; here, slow moving, deep biodiffusers, a situation which ultimately could threaten community resilience over the longer-term [1,34]. Further, this narrowing of trait expression may also generate a secondary loss of biodiversity resulting from changes in trophic interactions [108], facilitation [109], behaviour [25] or trait expression [12]. Indeed, the discrepancies between taxa richness and diversity and level of functioning call for a rethinking of management priorities that mainly focus on maintaining biodiversity levels [110,111].
Whilst innovative practical solutions to improve the condition and resilience of ecosystems remain in their infancy for marine benthic systems [112], interventions aimed at sustaining or improving ecosystem functioning are being applied without appreciation of the cascading effects these may have across multifunctional ecosystems [113,114]. Progression in this area will require better understanding of the mechanistic basis of trait expression and how species traits covary with the risk of extinction [9], whilst emphasising the importance of context and post-perturbation mechanisms in determining the functional consequences of altered biodiversity [115]. Rather than management or conservation efforts that primarily focus on the most biodiverse habitats, or the most common or productive species, our findings suggest that tailored management strategies that incorporate the dynamics of ecosystem responses to specific circumstances [99] and which enable the tracking of progress in recovery scenarios are likely to be of greater societal benefit [116].