Non-Linear Interactions Determine the Impact of Sea-Level Rise on Estuarine Benthic Biodiversity and Ecosystem Processes

Sea-level rise induced by climate change may have significant impacts on the ecosystem functions and ecosystem services provided by intertidal sediment ecosystems. Accelerated sea-level rise is expected to lead to steeper beach slopes, coarser particle sizes and increased wave exposure, with consequent impacts on intertidal ecosystems. We examined the relationships between abundance, biomass, and community metabolism of benthic fauna with beach slope, particle size and exposure, using samples across a range of conditions from three different locations in the UK, to determine the significance of sediment particle size beach slope and wave exposure in affecting benthic fauna and ecosystem function in different ecological contexts. Our results show that abundance, biomass and oxygen consumption of intertidal macrofauna and meiofauna are affected significantly by interactions among sediment particle size, beach slope and wave exposure. For macrofauna on less sloping beaches, the effect of these physical constraints is mediated by the local context, although for meiofauna and for macrofauna on intermediate and steeper beaches, the effects of physical constraints dominate. Steeper beach slopes, coarser particle sizes and increased wave exposure generally result in decreases in abundance, biomass and oxygen consumption, but these relationships are complex and non-linear. Sea-level rise is likely to lead to changes in ecosystem structure with generally negative impacts on ecosystem functions and ecosystem services. However, the impacts of sea-level rise will also be affected by local ecological context, especially for less sloping beaches.


Introduction
Climate change is expected to have significant effects on the marine environment through temperature change, ocean acidification, and accelerated sea-level rise, due to thermal expansion and melting of ice sheets [1]. Warmer sea temperatures are likely to cause changes in abundance, diversity and size composition of zooplankton [2,3], increases in abundance of southern invertebrate species and decreases of northern species in the northern hemisphere [4]. Acidification may have both direct and indirect impacts, via changes in the phytoplankton community, on bacteria and zooplankton [5,6,7]. Some of the greatest impacts of sea-level rise are likely to occur in intertidal sediment habitats. Intertidal habitats will rarely be allowed to transgress inland due to the high value of real estate and land behind the existing high tide zone, so they will be squeezed between rising sea-level and hard coastal defences. Coastal squeeze leads to the loss of intertidal area, sedimentary shifts towards coarser particles and more reflective morphodynamic states, greater tidal velocity, changes in water depth and (for estuaries) salinity, and increased storm surges [1,8,9,10].
Estuarine soft sediment habitats and their associated biodiversity are of great functional importance for the entire marine ecosystem and for human wellbeing in terms of the ecosystem services they produce, including high primary productivity, high secondary productivity, nutrient cycling, climate regulation, pollution control, decomposition, biodegradation and recreation [8,11,12,13,14,15,16,17]. They also provide nursery grounds for marine fish, and feeding and breeding areas for migratory birds and other species [18]. Many of these services are underpinned by ecosystem processes occurring within the benthos [19,20,21]. Understanding the likely impacts of the physical constraints imposed by sea level rise on benthic assemblages is therefore important for evaluating and managing threats to future ecosystem service provision.
Coarser sediment size, steeper beach slopes and higher exposure to wave action are predicted for intertidal areas of the typical Vshaped estuaries that characterise much of the European coastline [9,22]. Changes in these physical beach characteristics are likely to have an impact on benthic invertebrate biomass and body size distribution [22,23,24], and consequently on ecosystem processes. However, the nature of these impacts is likely to be determined by the interactions among the physical constraints of particle size, beach slope and exposure, collectively known as beach morphodynamic state [23,25], rather than any of these factors acting alone. At one end of the morphodynamic range, dissipative beaches with finer particles and gentle slopes have been shown to support a higher number of species and a more abundant macrofauna [26,27,28,29], whilst the opposite is true for reflective beaches at the other end of the range. Meiofauna are likely to be less affected than macrofauna by changes in exposure and sediment particle size [30].
Here we extend previous considerations of the potential impacts on estuarine benthic organisms of physical changes associated with sea-level rise to the ecosystem processes with which they are involved. We examine the relationships between abundance, biomass, and community metabolism of benthic organisms with beach slope, particle size and exposure, using samples that span a range of conditions from three very different locations in the UK. We test the following specific hypotheses [26,27,28]: (1) Total invertebrate abundance, biomass, and species richness will decline in response to increases in average sediment grain, steeper beach profiles and greater exposure to wave action; and (2) these changes will also be reflected in changes in community metabolism, a measure of ecosystem functioning. We also specifically consider the interactions between the different physical constraints and local ecological context in determining the overall impacts on beach fauna.

Ethics Statement
No specific permits were required for the field studies. All the sites sampled are public beaches, except for the Ythan which is a National Nature Reserve. Samples were taken from the Ythan under the permit issued by Scottish Natural Heritage to OceanLab, University of Aberdeen (Dr Martin Solan). None of the field studies involved endangered or protected species.

Study Sites and Sampling Method
We chose three estuaries within the UK: the Humber (from 53u339200 to 53u349500 N, and from 0u009400 to 0u039200 W), the Ythan (from 57u189500 to 57u209200 N, and from 1u599200 to 2u019 100 W) and the Firth of Forth/the Forth estuary (from 55u579300 to 56u009300 N, and from 3u069300 to 3u319000 W) on the basis of their wide geographical spread within the UK and previous history of research. We selected five or six sampling stations on each of these estuaries to provide a range of particle sizes, slopes and exposures (Fig. 1). The sampling stations were restricted to a defined section of the outer (marine) salinity gradient (.30, spring high tide) to minimise any potentially confounding effects of salinity [31]. At the mid-tide level at each sampling station, to correspond with maximum invertebrate abundance, a cylindrical core (10.4 cm diameter) was pushed into the sediment to the depth of 10 cm on a randomly chosen surface to sample macrofauna, meiofauna and sediment, thus ensuring that all variables were measured at a similar scale (grain). Sampling stations were established using a map and previous surveys [31,32,33] Four widely-spaced replicates were taken for macrofauna, one sample for meiofauna and one sample for sediment. No replicates were taken for meiofauna as the sampling unit (86 cm 2 ) is far greater than the scale of patchiness for this group (c. 3-10 cm 2 ) [34], such that variance between samples of this size is extremely low.

Biological Measurements
Macrofauna were separated from sediment using a 500 mm mesh, preserved in 70% ethanol, identified to species level wherever possible, and counted, using a low-power microscope. Meiofauna were separated from sediment using a 64 mm mesh, preserved in 70% ethanol and stained with Rose Bengal (5 mg l 21 70% ethanol). Because of the large core size used for sampling meiofauna and the large number of individuals collected, meiofauna from each site were sub-sampled after extraction and homogenisation. Meiofauna in these sub-samples ( = 1/200 th of core) for each site were identified to the lowest possible taxon, counted and measured to provide body dimensions for the calculation of body mass (Table S1).
Calculation of body mass. Body dimensions of individual animals were measured under a high-power microscope, and converted to body mass (dry weight) using established relationships ( [35]; Table S1). Where densities were high, body dimensions were measured and body mass calculated for a sub-sample, and then scaled up to the number in the full sample.
Community metabolism. Community metabolism was calculated as oxygen consumption by all individual invertebrates using data on mass-specific oxygen consumption factors derived by Banse [36], as used by Gerlach, Hahn, et al. [37]. Abundances of macrofauna and meiofauna were expressed as numbers core 21 , biomass as Ash Free Dry Weight (AFDW) as mg core 21 , and oxygen consumption as mlO 2 hr 21 core 21 . Thus, Log 10 (Oxygen Consumption) = 20.2599*Log 10 (AFDW) + 0.51 for macrofauna; Log 10 (Oxygen Consumption) = 20.24*Log 10 (AFDW) + 0.0096 for meiofauna; and Log 10 (Oxygen Consumption) = 20.2605*Log 10 (AFDW) 2 0.4424 for foraminifera [34,35]. Sea-Level Rise Impacts on Estuarine Processes PLOS ONE | www.plosone.org Beach Physical Characteristics Sediment particle size. Sediment cores were oven-dried at 80uC for at least 24 hours, after which the sample was thoroughly homogenised. Particle size composition was determined by dry sieving the material through a tower of sieves of mesh 4.75 mm, 2.8 mm, 1 mm, 500 mm, 250 mm, 125 mm and 64 mm. Median particle size (mm), 25% and 75% quartiles (Q 25 , Q 75 in mm) and the sorting coefficient QDQ~Q 25 {Q75 2 were determined by standard probit analysis. Silt content was taken as the weight of material passing through a 64 mm mesh by wet sieving.
Beach slope. The slope at each sampling station was calculated as: where H = vertical height (m) difference between the low and high shore, and D = distance (m) down the beach between the low and high points. Height difference was measured using a dumpy level and a staff, and distance using a surveyor's tape. Because the values for slope on these beaches are very small, they are presented in the text as slope610 3 , to facilitate inter-site comparisons.
Exposure. Exposure was calculated using a modified form of the Thomas Exposure Index (TEI) [38], derived from wind velocity, direction, duration and the effective fetch: where W = (percentage of time wind blows in a sector (22.5 degree sectors of the compass rose)/100) 6 (mean wind speed (kn) 2 ). F = Fetch in nautical miles (100 nautical miles maximum). CS = Extent in nautical miles of water ,6 m deep adjoining the shore. DS = Extent in nautical miles of water within the fetch ,6 m deep, but not adjoining the shore.
W was calculated using the wind data from windfinder.com, and F, CS, and DS were derived from UK Admiralty charts.

Statistical Analysis
Statistical models (mixed model for macrofauna and linear regression with GLS extension for meiofauna) were developed [39,40] for abundance, biomass and oxygen consumption of macrofauna and meiofauna, in relation to the physical beach variables. Physical beach variables (median particle size, beach slope and wave exposure) were measured for each station, and four replicates were included in the model for macrofaunal abundance, biomass and oxygen consumption and no replicates were included for meiofaunal abundance, biomass and oxygen consumption. Silt content was not included in the models as it was strongly correlated with median particle size, and the sorting coefficient QDQ was not significant in any models. We included the three estuaries initially as a random local ecological effect in all models; this term was highly significant in the model for macrofauna, although not in that for meiofauna. For the macrofauna model, we therefore used a mixed modelling approach, including median particle size, beach slope and exposure as physical fixed effects and estuary as a random local effect. For the meiofauna model, we included median particle size, beach slope and exposure as effects, and used a linear regression model with a GLS extension to allow for unequal variance associated with median particle size [41,42]. All the possible two and three way interactions were included in both macrofaunal and meiofaunal models.

Beach Physical and Biological Characteristics
The Humber and the Forth estuaries had a similar range of particle size (with the Ythan range slightly lower), and the Humber had much higher exposure than the Ythan and the Forth (no overlap in range). The Ythan and the Forth also had a similar range of beach slopes, whereas those for the Humber were much shallower ( Table 1, Fig. 1).
The Ythan and the Humber had similar numerically dominant species and species richness of macrofauna, but differed markedly in the numbers of individuals recorded ( Table 2). The Forth had intermediate abundances, but a greater range of species, and was dominated numerically by polychaetes. The Forth had the most meiofaunal taxa represented, including mites (Acarina) and archiannelids. Nematodes and foraminiferans (live as opposed to dead shells) occurred in large numbers, especially on the Humber and the Ythan, and Copepods and Turbellarians were abundant on the Ythan.

Relationships of Macro-and Meiofauna with Beach Physical Characteristics
Macrofaunal abundance, biomass and oxygen consumption. Median particle size and beach slope had a similarly strong influence on the abundance of macrofauna (Table 3), whilst exposure had relatively less influence on abundance. The three-way interaction term particle size 6 beach slope 6 exposure was also significant. Sediment particle size, beach slope and wave exposure are usually associated to each other (the beach with high wave exposure tend to have steeper beach profile and coarser sediment). Because the local ecological effect (estuary) was highly significant, the graphs to show the predictions of the mixed model for macrofaunal abundance were made for minimum, mean and maximum exposure and beach slope for each estuary separately (Fig. 2). For shallow-sloping beaches, there were no consistent relationships between abundance and exposure or median particle size across the estuaries, and the local ecological effect dominated the patterns observed (Fig. 2). Thus, for shallow-sloping beaches, there was a weak negative relationship between abundance and particle size for the Humber stations, a slightly stronger, positive one for all the Ythan stations and no relationship for the Forth stations. However, for intermediate and steep beaches, there were clear and consistent non-linear relationships between abundance and exposure and particle size, which overrode the local ecological effect (Fig. 2). For the most sheltered beaches, there were strong negative relationships between abundance and particle size. However, as exposure increased, this relationship became more weakly negative and then positive at intermediate exposures, and then became negligible or very weakly positive at the highest exposures (Fig. 2). Macrofaunal biomass and oxygen consumption exhibited the same non-linear relationships with physical factors as macrofaunal abundance (Figs. S1 and S2). The most significant single factor affecting both biomass and oxygen consumption was beach slope (Table 3), followed by median particle size and exposure. The three-way interaction term particle size 6 beach slope 6 exposure was also significant for both macrofaunal biomass and oxygen consumption (Table 3). Biomass and oxygen consumption show contrasting trends because weight-specific oxygen consumption increases as organisms become smaller [36,37].

Meiofaunal
abundance, biomass and oxygen consumption. Beach slope had the greatest effect on meiofaunal abundance (Table 3), followed by median particle size and exposure. The three-way interaction term particle size 6 beach slope 6 exposure also had a significant effect on meiofaunal abundance; meiofaunal abundance decreased with increasing exposure and increasing particle size on the shallow-and  intermediate-slope beaches (Fig. 3). However, on the steep beaches, there was evidence of a non-linear response to the interaction between exposure and particle size, similar to that observed for macrofauna. Meiofaunal biomass and oxygen consumption exhibited the same non-linear relationships with physical factors as meiofaunal abundance (Figs. S3 and S4). Meiofaunal biomass and oxygen consumption were affected most significantly by exposure, followed by beach slope and median particle size ( Table 3). The three-way interaction term particle size 6 beach slope 6 exposure was a significant factor determining both meiofaunal biomass and oxygen consumption. The two-way interactions were also significant for all the models, but these lower term interactions are not discussed here.

Discussion
Our analyses revealed strong relationships between macrofaunal and meiofaunal abundance, biomass and oxygen consumption and beach slope, particle size and exposure. That relationships exist is not unexpected, but here we have been able to tease out the relative importance of the different variables and, more importantly, their linear or non-linear nature. The most significant factor affecting macrofauna was slope, followed by particle size and exposure. For meiofauna, the most significant factor was exposure, followed by slope and particle size. For the shallowest beaches, the local context (which estuary it was) overrode these other relationships, but for intermediate and steeper beaches, the physical factors were the dominant influences. The overall pattern for the relationship for intermediate and steep beaches was negative, with invertebrate abundance, biomass and oxygen consumption declining with increasing slope, particle size and exposure, but the effects were complex and non-linear, with the relationship switching at intermediate levels of exposure.
Previous studies have shown that dissipative beaches with fine particles and gentle slopes support higher numbers of species, in greater abundance and biomass [15,23,26,27,28]. The present study confirms these trends for the meiofauna, showing declines in abundance, biomass, and oxygen consumption with increasing particle size, beach slope, and exposure. However, our results go much further than previous analyses in demonstrating that under certain conditions these relationships are likely to be more  Table 3. L-ratios for regression models of abundance, biomass and oxygen consumption for macrofauna and meiofauna with respect to the single factors median particle size, slope and exposure (each df = 4) and the three-way interaction between them (df = 1). complex and non-linear. This has important implications for the impact of sea-level rise on ecosystem service provision such as reduced food resources, reduced water quality, reduced primary productivity and change in release of nutrient to the overlying water column, since the effects of steeper slopes, increasing particle size and increased exposure on ecosystem structure and function will be moderated by local conditions, especially for less sloping beaches.
McLachlan et al. [30] reported greater numbers of meiofaunal harpacticoid copepods and fewer nematodes in relation to larger median particle sizes. Those authors found that nematodes tended to dominate in finer sands (,330 mm), with overall biomass higher, similar to the present study, where all sites had particle sizes ,330 mm and were also dominated by nematodes. In contrast, Rodriguez et al. [43] found increased meiofaunal biomass with increasing wave exposure, but their sites were much more exposed. Taken together, these results emphasise the complexity of the relationships between biodiversity elements and physical constraints on beaches, which probably contribute to the non-linear trends observed by others for ecosystem service provision in coastal ecosystems [14].

Interactions between Beach Physical Factors and Local Ecological Context
McLachlan [25] and Raffaelli and Hawkins [23] argued that the interaction between particle size, beach slope and exposure is probably as important, if not more so, for predicting changes in beach fauna than the effect of any single physical factor. This is formally confirmed in the present study for the first time: the simple bivariate relationships depicted in the literature such as graphs of abundance or species richness against sediment particle size (e.g. [23], p59, Fig. 2.10) are clearly inadequate descriptions of the true, more complex relationships, and would be misleading if used to predict impacts of sea-level rise.
The meiofaunal models revealed linear relationships between variables, but all the macrofaunal models indicate non-linear relationships. Since maximum faunal abundance and diversity may change its exact location along the intertidal gradient depending on exposure, and the sampling stations were fixed at a geographical mid-tide level in this study, it could be argued that the observed non-linear responses of macrofauna with exposureslope could be due to sampling at this fixed point. However, we think this is unlikely given that the meiofaunal relations might then have been expected to be similarly non-linear. The underlying causes of the non-linear macrofaunal relationships thus remain unclear.
The location of the three estuaries was highly significant for models of abundance, biomass and oxygen consumption of macrofauna, but not for the meiofauna. This suggests that the importance of physical factors outweighs that of local context for meiofauna, although for macrofauna, this is only the case under the limited conditions of less sloping beaches. The exact nature of these local effects is unknown, but they must reflect other local variables which we did not measure, such as differences in predation by birds and fish, or differences in water quality.

Impact of Sea-level Rise on Ecosystem Services from Estuarine Systems
Accelerated sea-level rise is expected to make beaches coarser and steeper, and more reflective in morphodynamic state [10,29]. The predicted increase in the frequency of storm surges [1] will add to this effect. Our results suggest that if beach physical factors change as predicted, then sea-level rise will make estuarine intertidal areas less diverse and less productive, through declines in abundance, biomass, and community metabolism as well as through the loss of area due to coastal squeeze. Mechanistic understanding of the relationships between biodiversity, ecosystem functions and ecosystem services remains uncertain for estuaries [42], although experimental work on benthic biodiversity and ecosystem functioning suggests that reductions in macro-and meiofaunal biomass and metabolism are likely to have significant effects in reducing nutrient cycling within sediments and release of nutrients to the overlying water column [19,39,44,45,46]. Our results have highlighted the complex and non-linear nature of these relationships, and more precise prediction of the impact of sea-level rise on ecosystem service provision will depend on a much improved understanding of the interactive effects of local ecosystem contexts and globally-operating physical constraints on ecosystem structure and function. Figure S1 Predicted macrofaunal biomass based on the minimal adequate regression model for each estuary. The values on graphs show minimum, mean, maximum exposure values of each estuary. H, Y, and F stand for the Humber, the Ythan, and the Forth, respectively. Shallow, Intermediate, and