Ecosystem indicators: predicting population responses to combined climate and anthropogenic changes in shallow seas

(e


Introduction
In the well-studied North Sea, temperature rises and other physical changes like ocean acidification, reduced oxygen and sea-level rise are having an impact throughout the food web, with effects on abundance, distribution and biodiversity seen in plankton (Edwards et al. 2020), fish (Wright et al. 2020), seabirds (Mitchell et al. 2020) and mammals (Evans and Waggitt 2020).One of the more likely solutions to combat climate change is the introduction of large-scale offshore renewable energy (ORE) developments (wind, tidal and wave) of 100s of gigawatts (GWs) (IRENA 2019).However, the introduction of new structures and the extraction of more energy will have cumulative effects within the world's shallow seas and therefore will also influence whole ecosystems (Christiansen et al. 2022, Daewel et al. 2022, Dorrell et al. 2022).The size of these developments may also end up using more than 30% of shallow sea space, and some of the consequences that will most likely follow include physical habitat change (e.g. increase in stratification of the water column) (De Dominicis et al. 2018), displacement of fisheries (Kafas et al. 2017) and possible creation of de facto marine protected areas (MPAs) (Raoux et al. 2019).There is an urgent need for a better level of understanding of the combined effects of climate and anthropogenic change if we are to use marine spatial planning (MSP) and de facto MPAs effectively (Gissi et al. 2019).Physical processes controlling rates of mixing and stratification are of key importance to understanding changes in the level of primary production on a relatively fine to local scale (Schultze et al. 2020) but also on an ecosystem-wide spatial scale (Ludewig 2015).Thus, understanding how changes to stratification due to the introduction of ORE developments may change primary production and lead to further impacts on ecosystem services such as fisheries, but also on mobile marine predator populations, is essential.
However, to address that, we first need a much greater understanding of how different marine habitats and their associated ecosystems, and specifically their multiplicity of physical and biological interactions, are likely to change across different large-scale habitats in space and over time with both climate and anthropogenic transformations.To address this, we first need to identify good indicators of habitat and ecosystem change.That requires an assessment of what makes a good indicator of change in marine ecosystems and what exactly the indicator is telling us about the changes across space and over time; and, therefore, provide early warning for management decisions and 'best use' practice under future climate and increasing anthropogenic use scenarios.
One way to understand the drivers of ecosystem change is to use ecosystem models with time series of multispecies population characteristics, as well as both physical and biological ecosystem components, such that patterns of species population change can be quantified across space and over time, under different climate and/or anthropogenic scenarios.
Understanding the drivers of ecosystem change is challenging because of the variability in observations, for example due to imperfect methods of observation (Link et al. 2012) and uncertainty in potential associations due to external forces like climate change (Proulx et al. 2005).However, probabilistic methods such as Bayesian networks (BNs) can be used to capture ecological patterns between variables (Hui et al. 2022) and reveal trends across space and over time (Tucker and Duplisea 2012), without requiring specific information on mechanisms and huge amounts of observational data, in contrast to other ecosystem models (Uusitalo 2007).Such probabilistic models allow predictions to be made across different spatial and temporal scales and with a range of indicator species or functional groups representing all trophic levels.In particular, BNs use a network structure and inference that allow us to ask, 'what-if?' type questions of the data.For example, what is the probability of seeing a change in the biomass of pelagic fish assemblage, given that we have observed a change in the levels of fishing, as well as a change in the stratification of the water column caused by ORE developments, and a change in the bottom temperature from climate change?
This investigation uses dynamic hidden BN ecosystem models that were developed from identified consistent physical and biological indicators that were found to indicate patterns of ecosystem change and predict the ecosystem state within four regions with different habitat types (shelf edge, oceanic influences, deep and shallow seas) in UK waters over a 30-year period (Trifonova et al. 2021).An indicator is a variable that has been found to help predict/indicate patterns of habitat and/or ecosystem change over the last 30 years within the North Sea.We examine the models' accuracy in terms of their ability to reproduce observations of the trends (increases versus decreases) in all the ecosystem components (oceanographic processes as well as species/functional groups at all trophic levels).Therefore, we can explore in more detail whether the different regions show similarity in population trends given the dynamic aspect of the multiplicity of interactions across all trophic levels and their stressors at the level of the ecosystem.Then the ecosystem models were used in combination with 'what-if?' scenarios based on potential changes in climate (increase in temperature), climate and/or ORE developments (increase in stratification), and anthropogenic change (increase versus decrease in fishing) to explore the trends of either increases or decreases of different ecosystem components in response to such changes.By providing an understanding of the reactive responses across all trophic levels, we also aim to discover patterns of the true dynamic nature of bottom-up versus top-down effects across trophic levels and within regions with contrasting habitat types.Therefore, we can identify whether the same biophysical variables are found to be consistent indicators of ecosystem-level changes across space and over time.Here, by identifying indicator species and habitats, we further build confidence in our ability to examine potential impacts of bottom-up versus top-down effects that could inform strategies for coping with and adapting to climate change, the placement of ORE developments and potential knockon effects, such as fisheries displacements.Being able to predict implications of climate and other anthropogenic changes on ecosystem components and their relative adaptability and resilience will be useful to guide what species and/or habitats are more resilient/at-risk to what type of disturbances, and therefore what management decisions are required for the future sustainable management of marine ecosystems between different uses of our shallow shelf seas and under the influence of future changes.

Habitat type
The choice of regions was based on both physical and biological characteristics.Four spatial regions were defined: Shetland/Orkney, west of Scotland, and deep and shallow central North Sea regions (Fig. 1).We wanted to compare a range of different habitats (shelf edge, Atlantic influenced, deep and shallow seas) within the regions that are key features of interest with respect to future ORE developments, with some regions more likely to have large-scale development of either wind (floating in deeper (> 50 m), and static in shallow North Sea regions, < 50 m), wave (Atlantic influenced with the largest fetch) or tidal (Shetland/Orkney region).

Ecosystem components
The time series input data consisted of annual values  from the summer season (July, August, September and October) as either mean values of physical variables (e.g.temperature) or cumulative values of biological variables (e.g.net primary production), or maximum values of physical and biological variables: current speeds and maximum Chla, respectively.Biological variables for population dynamics included total annual abundance, biomass or mean breeding/pupping success per spatial region (Table 1).We refer to all the variables in the study as 'ecosystem components' but distinguish components based on them being either physical (e.g.horizontal currents speed) or biological (e.g.sandeel larvae) indicators.The ecosystem components in the study were chosen as they cover the main physical and biological variables that have been shown to be important to marine mammals and seabirds and their prey (Carroll et al. 2015, Wakefield et al. 2017, Chavez-Rosales et al. 2019).These will alter with climate change, and also with the next biggest change to our shallow seas: very large extraction of energy from ORE (van der Molen et al. 2014, Wakelin et al. 2015, Holt et al. 2016, Sadykova et al. 2017, Boon et al. 2018, De Dominicis et al. 2018).See Supporting information for more details on the choice of ecosystem components.

Model description and 'what-if?' scenarios
The modelling approach is a dynamic BN model with a hidden variable (HDBN) that is a modified version of the model developed in Trifonova et al. (2015Trifonova et al. ( , 2017)).See Supporting information for more information on modelling time series with dynamic BNs.Each model (four models per region) was built from the identified high dependency interactions for capturing similarities in the temporal trends of the paired ecosystem components to increase certainty of predictions of population changes.In this way, the region-specific indicators that were found to be the most consistent (i.e. the ecosystem components that were most often found in the paired high dependency interactions) for predicting the ecosystem state and indicating patterns of the associated dynamics for each spatial region and their changes over time were identified (Trifonova et al. 2021).Therefore, each HDBN ecosystem model captures the spatial and temporal variability of multiple biophysical interactions throughout the trophic chain, ensuring that the strongest relationships (i.e.relationships of high dependency that are predictive in an informative, not causal aspect), and so the most consistent indicators of ecosystem change, are the ones identified in this process.From the strongest relationships, up to three indicators (i.e.'parent' nodes) were selected that drive the target ecosystem component (i.e.'child' node) and were used to build the region-specific modelling structures.For example, in Fig. 2b, sandeel larvae (the 'child' node) is driven by two 'parent' nodes: BT and fisheries catch.
The model incorporates lower tropic levels (e.g.Chl-a, NetPP, zooplankton assemblages), commercial and important prey fish groups (sandeel larvae, pelagic and demersal) and top predators (seabirds and marine mammals).The model also incorporates physical indicators (BT, SST, HSpeed, VSpeed and PEA, Table 1) which are those that will produce large bottom-up changes with increasing climate change and the large-scale ORE developments through temperature, mixing and stratification increases (Sharples et al. 2020, Dorrell et al. 2022).The model is run in the four contrasting habitat types (Fig. 1) where the placement of new structures and effects of energy extraction of ORE developments can also change temperature and levels of stratification (Christiansen et al. 2022) and where varying levels of fishing will be used to explore expected top-down effects on the ecosystem components (Lynam et al. 2017).
Hence, we perform 'what-if?' scenario analyses, following moderate and extreme changes to mean bottom temperature (BT) or mean sea surface temperature (SST) (increase by 1°C and increase by 2°C) and mean potential energy anomaly (PEA) (increase by 50 J m -3 and increase by 100 J m -3 ), under three levels of fishing: high, medium and low to investigate the effects on the ecosystem components.Then, the ecosystem components were classified as either consistent indicators of bottom-up effects driven by physical changes in BT/SST, and/or PEA versus indicators of top-down effects driven by fishing or indicators of both effects across space and over time.This was determined based on whether any type of change (increase or decrease) was consistently predicted for an ecosystem component in response to the scenario changes in BT/SST, PEA (i.e.bottom-up driven), and/or fishing (i.e.top-down driven) across space and over time.To identify whether an indicator is bottom-up or top-town driven, we also used the responses of other ecosystem components to 'what-if?' scenarios that drive the population trends.In this way, based on how often the type of indicators (bottom-up versus top-down or both) were identified across regions, we were able to define if regions were more or less often driven by the same effects.In addition, to identify whether a region is bottom-up, top-down driven or both, we considered the region-specific indicators that were found to be the most consistent (i.e. the ecosystem components that were most often found in the paired high dependency interactions) for building the models.
Non-parametric bootstrap (re-sampling with replacement from the training set, Friedman et al. 1999) was applied 250 times to each baseline model and scenario run to obtain validation in the predictions for each region.Model performance, in terms of sum of squared error (SSE), was assessed for each baseline model and predictions were compared on a year-to-year basis versus the original input data: The modelling structure varies depending on the spatial region (e.g.Shetland/Orkney model shown in Fig. 2b), but the general form is presented in Fig. 2a.Each model included a single hidden variable that was modelled as a discrete node with two states.Given the probability distribution over X[t] where X = X 1 …X n are the n variables observed along time t, to predict each ecosystem component, we inferred the component at time t by using the observed evidence (or available data) from t − 1.We used an exact inference method: the junction tree algorithm (Murphy 2002).The data were standardised prior to conducting the experiments to a mean of 0 and standard deviation of 1.
See Supporting information for more information on the scenarios and experiments.

Model performance and accuracy: comparison of sum of squared error and timeseries of population trends
The dynamic Bayesian network model with a hidden variable (HDBN) across all the regions reported outcomes of high predictive accuracy for the ecosystem components by using  up to three indicators (see symbols in Table 2) that generally provided low values (less than 25.00) of SSE (68% of the ecosystem components having values lower than the threshold of 25.00, Table 2).We found the threshold of 25.00 to be most appropriate based on examining the range of SSE values across habitats as well as across ecosystem components within a habitat.Comparison of the predictive performance of the region-specific ecosystem models indicates varying spatially predictive accuracy.The Shetland/Orkney region had the most ecosystem components predicted with values of SSE less than 25.00 (86% of components), followed by the two regions within the central North Sea (66% of components) and finally, west of Scotland region (53% of components).
Within the majority of ecosystem components that were accurately predicted, sandeel larvae, NetPP and the two seal species were predicted with high accuracy across all regions, whilst the A4 zooplankton, pelagic fish assemblages and kittiwake were predicted most accurately in three out of the four regions.However, there were 32% of the components with lower accuracy.For example, the A5 zooplankton assemblage in the deep central North Sea had the highest SSE of 29.71 (Supporting information).We show the time series of the predicted ecosystem components from the baseline HDBN models to visually demonstrate how well they performed in reproducing the inter-annual variability and long-term patterns (always shown as blue lines) versus the original input data (black lines) (Fig. 3-6).Note, we only show some illustrative examples, with the time-series for all remaining ecosystem components and their predicted responses to the scenarios, including 95% confidence intervals calculated from the bootstrap predictions' mean and standard deviation, shown in the Supporting information.The HDBN models were able to capture many of the changes (increases versus decreases) of the ecosystem components across space and over time, predicting the general trends in population dynamics for all functional groups and species using three or fewer indicators.In the figures (Fig. 3-6), the solid black symbols represent the indicators identified as the most confident data-driven interactions (Trifonova et al. 2021), and were therefore used  to build the region-specific HDBN modelling structures.We encountered time lags in some of the predictions (e.g.guillemot in Fig. 6c), which we believe is to be explained by structural uncertainties (i.e.species-specific relationships and/or colony-specific drivers used to build the models, see Discussion on bottom-up and top-down indicators).Although we encountered some time lags, the models performed well in reproducing the inter-annual variability and long-term patterns across space and over time (e.g.guillemot and kittiwake annual breeding success declines during late 1990s-early 2005 and increase after 2015).

'What-if?' scenarios
Next, we describe the results from the scenarios by examining if the predictions, and specifically if the temporal trends of the ecosystem components, were predicted to increase or decrease in comparison to the original input data over the same time period, separately for the different regions.The scenario analyses included changes in the indicators: bottom temperature (BT), sea surface temperature (SST), potential energy anomaly (PEA) and levels of fishing (denoted by coloured symbols in Fig. 3-6).Across all regions, however, whether an increasing or declining trend would be predicted for an ecosystem component, following the scenarios, was dependent on the component-specific interactions with the physical, biological and/or fishing indicators used to build the HDBN models in the different regions (denoted by solid black symbols).That is why, in some figures, coloured symbols indicating scenario changes are either temperature, PEA and/ or both dependent on the component-specific interactions in that region.
We summarize the predicted trends (increasing or decreasing blue coloured arrows) for all ecosystem components in response to the 'what-if?' scenarios in Table 3.To determine the illustrative examples in Fig. 3-6, we considered the values of SSEs and timeseries plots, thus showing the top one and/or top two ecosystem components with the lowest SSE values per region.Also, the selection of ecosystem components, shown in the figures below (Fig. 3-6), was based on whether the ecosystem components were identified as consistent indicators of bottom-up effects driven by physical changes in BT/SST, and/or PEA versus top-down effects driven by fishing or indicators of both effects across space and over time.

Bottom-up indicators: lower trophic levelsmaximum Chl-a, net primary production and zooplankton
Following scenario increase in BT and/or PEA, both Chl-a and NetPP were predicted to have a declining trend over time in the deep central North Sea, and Shetland/Orkney, respectively (Fig. 3b,d).In the deep central North Sea region, the majority of zooplankton assemblages were predicted to have a declining trend, following scenario increases in BT and/PEA.For example, a declining trend in the deep central North Sea region was predicted for A4 (Fig. 3f ), following scenario increases in PEA, which also predicted declining trends for Chl-a (Fig. 3b) and sandeel larvae (Fig. 5b) that drive the A4 annual abundance trend in this region (Fig. 3e).

Bottom-up indicators: top predators -gannet, harbour porpoise and harbour seal
The upper trophic level species that were found to be predicted with either increasing or decreasing trend due to scenario increases in BT or PEA were gannet, harbour porpoise and harbour seal.Gannet breeding success in the west of Scotland was predicted to be increasing over time due to scenario increases in PEA only (Fig. 4b).Note, there were no explicit changes of either increase or decrease in the other ecosystem components that were shown to drive gannet annual breeding success (Fig. 4a), such as A5 and sandeel larvae abundance, which is, in turn dependent on fishing levels (Supporting information).Harbour porpoise in the Shetland/Orkney region, were not predicted with an explicit trend of either increase or decrease (Fig. 4d) given  the scenario increases in BT and PEA.For the zooplankton assemblages (A4, A5 and A6) that drive harbour porpoise annual abundance (Fig. 4c) there were no explicit trends of either increase or decrease for A5 but there were increasing trends for A4 and A6 due to scenario increases in BT and PEA (Supporting information).The trend for harbour seal abundance was predicted as a declining trend (Fig. 4f ), due to scenario increases in BT in the west of Scotland region.Harbour seal annual abundance in that region is also driven by Chl-a and NetPP (Fig. 4e), with the scenario predicting declining trends (Supporting information).

Top-down indicators: sandeel larvae and fish assemblages
The species/groups that were found to be predicted with either increasing or decreasing trend due to scenario changes in fishing were, unsurprisingly, sandeel larvae and commercial fish species.Sandeel larvae abundance was predicted to decline, following a scenario increase in BT and a high fishing level in the deep central North Sea (Fig. 5b).The pelagic and demersal fish assemblages were predicted with a declining trend following the scenario of high fishing levels (Fig. 5d,f ) in the shallow central North Sea.Following scenarios of low and medium fishing levels, an increasing trend was predicted, that was more explicit for the demersal fish assemblage (Fig. 5d), in comparison to the pelagic fish assemblage (Fig. 5f ).For sandeel larvae that drive fish assemblages' annual biomass trends (Fig. 5c,e) in that region, varying trends were predicted that were dependent on fishing levels (Supporting information).A decreasing trend was predicted for the zooplankton assemblage A2 that drives the pelagic assemblage (Supporting information), whilst no explicit trend of either increase or decrease was predicted for NetPP (Supporting information) that drives the demersal fish assemblage annual biomass trend (Fig. 5c).

Bottom-up and top-down indicators: top predators -grey seal, guillemot and kittiwake
There were some upper trophic level species (grey seal, guillemot and kittiwake) that were found to be predicted with either increasing or decreasing trend due to scenario changes in BT/SST, PEA and fishing, thus identifying those species Table 3. Predicted trends: increasing or decreasing (denoted by directed arrows) for the ecosystem components per region.Double directed arrows indicate both increasing and decreasing years rather than an explicit trend of one direction.The symbols denote component-specific interactions that are used to build the Bayesian network with hidden variable (HDBN) models across the four regions.PEL, pelagic fish assemblage; DEM, demersal fish assemblage.Page 13 of 18 as indicators of both bottom-up and top-down driven effects across regions.Grey seal productivity in the deep central North Sea was predicted with an increasing trend until 2010, but after that, the trend was predicted to decline (Fig. 6b).The predicted trend was dependent on the scenario increases in BT but also on the levels of fishing that led to changes in the predicted trends for sandeel larvae abundance (Fig. 5b) that drive grey seal annual productivity in that region (Fig. 6a).Guillemot breeding success in the Shetland/Orkney area was predicted to have a declining trend in the late 1990s-early 2000s, compared to the original input data; however, after 2005, the trend was not predicted with an explicit increase or decrease, nor did the predicted trend match any of the decreases (e.g. in 2008, 2011) of the original input data (Fig. 6d).The predicted trend was dependent on the scenario increases in BT and PEA, which predicted increasing trends for the zooplankton assemblages A2 and A6 (Supporting information), and fishing levels, which predicted varying trends for the demersal fish assemblage (Supporting information) that drive guillemot annual breeding success in that region (Fig. 6c).Kittiwake breeding success in the west of Scotland was not predicted with an explicit trend of either increase or decrease (Fig. 6f ).Similarly to guillemot, declining trends were predicted during the late 1990s-early 2000s; and then, after 2015, the predicted trend did not match any of the decreases during 2005-2008, compared to the original input data.The predicted trend was dependent on the scenario increase in SST and fishing levels, which predicted varying trends for sandeel larvae (Supporting information) that drive kittiwake annual breeding success in that region (Fig. 6e).

Model performance and accuracy: habitat matters
Our results highlight the need to include region-specific ecosystem level changes and dynamics of the multiplicity of interactions when building predictive models of complex and heavily exploited ecosystems within shallow seas, such as the North Sea.Hidden variables can pick up ecological patterns in the data that agree with ecosystem change that might not be strictly represented within the model structure (Uusitalo et al. 2018).Indeed, the accurate performance of the models was likely due, at least in part, to the inclusion of the hidden variable that reduced the likelihood of introducing spurious interactions into the analysis and allowed for more plausible modelling network structures.Overall, we found variability in species response to change driven by the habitat type (shelf edge, oceanic influences, deep and shallow seas) associated within the different regions.Therefore, effort should be emphasized on the integral habitat and, specifically, to determine what aspects of the habitat (e.g.bottomup forcing through physical trends and primary production dynamics versus top-down human pressures) are relevant to marine animals to be able to address a wider range of ecological questions.With the improved understanding of the exact bottom-up (e.g.levels of mixing and stratification) versus top-down (e.g.predators and fishing) mechanisms that influence habitat use by marine animals across spatial (< 1 km through to 1000 km) and temporal (days through to years) scales, the effects of biophysical interactions on populations and ecosystems and how these vary with climate change can be better understood.Therefore, to understand and predict ecosystem level changes caused by large-scale ORE developments, as well as separating out the effects of climate change, one needs to start with an appreciation of the changes that both these factors have on bottom-up forcing of levels of mixing (and stratification), and primary production, an overall indicator of ocean health (Tett et al. 2008(Tett et al. , 2013)).This new level of understanding will require the more targeted collection of at-sea data with simultaneous information on both animal behaviour and physical processes at the appropriate physical scales.Our findings of some less accurate predictions (e.g. the A5 zooplankton assemblage in the deep central North Sea, Table 2) have indicated the need to investigate the potential for different units of measure (mean, median or percentile) for some physical and biological indicators.Although the mean is the most widely used unit across studies, examples such as maximum Chl-a (Scott et al. 2010) illustrate that other units can have a more meaningful interpretation, when the specifics are understood of why the variable is important to mobile marine species.As such, it is imperative that this type of study is performed in habitats that have spatially and temporally rich baseline data, such that the understanding of mechanisms is achieved and can be subsequently transferred to other locations and across larger regional and shelf-wide scales through the use of ecosystem modelling.

Regional summary: bottom-up dominated regions more robust than top-down regions
Across the four spatial regions, there were a lot of similarities in the trends (increasing versus declining) of species response to change.However, it was noted that the species in the deep and shallow regions within the central North Sea were predicted more often by declining trends, whilst the Shetland/Orkney region was the only region in which many species were predicted with increasing trends, followed by the west of Scotland region.These results suggest that the ecosystem dynamics in the latter two regions might be more robust to natural and/or anthropogenic changes.To further support the potential higher level of robustness in the latter two regions, previous results on the modelled ecosystem state by the hidden variable also suggested that their dynamics were in a more desirable (predictable) state, with the timing of the North Sea 'regime shift' more clearly identified for the central North Sea regions, rather than for the other two regions (Trifonova et al. 2021).The Shetland/ Orkney region is strongly dominated by bottom-up processes driven by BT and PEA, but with changes in fishing levels also found to have strong top-down effects.This region has large resources of tidal power; however, the extraction of reasonable amounts of tidal energy (6 GW) has been shown to have effects on PEA tens to hundreds of miles downstream (De Dominicis et al. 2018).For the other oceanic influenced region, the west of Scotland, which is the type of region targeted for wave and floating wind, our results suggested that bottom-up driven effects on the dynamics are also dominant.These findings suggest the possible mechanisms behind region-specific and temporal ecosystem changes are needed to reveal insights into the details of marine ecosystem structure and function under the stress of both climate change and anthropogenic forces.
Stronger bottom-up effects driven by BT and PEA also dominated the deep central North Sea region; however, topdown effects driven by fishing (captured well in sandeel larvae abundance predictions) were also found to be highly important for the ecosystem dynamics.The number of interactions that direct the model for this region were dominated by more biological indicators (e.g.sandeel larvae, Chl-a, Table 2), in comparison to all the other regions.These results suggest that both bottom-up and top-down forcing can alter the ecosystem dynamics for this region through species interactions mediated by key species (sandeel larvae) and/or processes (primary production), and that it is the strong variation in the relative importance of bottom-up versus top-down forcing over time that characterizes this region.
Conversely, the shallow central North Sea region was characterised by stronger top-down effects, with only some bottom-up effects, driven by SST and PEA.Predicted declining trends dominated the region's ecosystem dynamics, which we suggest is most likely due to large influences from densitydriven stratification, high turbidity and rapid warming (van Leeuwen et al. 2016), in comparison to the other regions.This characterized the region as more temporally variable (Capuzzo et al. 2018) and an 'early indicator' of profound ecosystem changes that may later be seen within the northern and deeper regions, given future climate and anthropogenic changes.By investigating baselines, our findings indicate that we need to pay much closer attention to how different habitat types of our global shallow seas are managed, as the same human (large-scale ORE developments, changes in fishing, etc.) and climatic changes will be felt very differently in different regions.

Bottom-up indicators: lower trophic levels -maximum Chl-a, net primary production and zooplankton
The biological ecosystem components (Chl-a, NetPP and zooplankton assemblages) were identified as indicators of bottom-up effects, as expected, as plankton abundance is known to be directly driven by the physical forcing of temperature and stratification (Simpson and Sharples 2012).In general, where Chl-a and NetPP were predicted with decreasing regional trends across all scenarios, the population trends of higher trophic levels were also predicted to decline, suggesting that ecosystems will respond consequently to anthropogenic changes in both temperature and stratification.Notably, NetPP was more often predicted to be driven by vertical current speed with field data also indicating that mean seasonal change in NetPP is through vertical mixing (Zhao et al. 2019).This is an important distinction, as NetPP values are more highly influenced by depth, with shallower, more mixed areas generally having higher depth integrated production; whilst Chl-a, which has a critical habitat role in shelf seas (Mérillet et al. 2020, Trifonova et al. 2021), can be influenced by a greater variety of weather and physical variables (Molinero et al. 2013).The potential for spatial changes in NetPP can have strong implications for fisheries production given future climate and fisheries management changes (Capuzzo et al. 2018).Therefore, understanding the drivers of primary production changes in response to ORE developments across regions of contrasting depth and hydrodynamic conditions is important for successful ecosystem-based management, and could influence the implementation of MPAs and highly protected marine areas (HPMAs) (Benyon et al. 2020).
The modelled outputs showed that predicted changes in zooplankton abundances consequently progress to affect population trends of fish assemblages (Fig. 5) and even top predators (Fig. 6), as it has been shown from direct analysis with abundance changes in zooplankton assemblages affecting fish recruitment (Beaugrand et al. 2003).Predicted zooplankton assemblages showed mostly decreasing trends within the North Sea regions but with increasing trends predicted for A2 (e.g.Calanus helgolandicus) and A6 (e.g.Calanus finmarchicus) in the Shetland/Orkney region (Supporting information).These findings again confirm the importance of understanding regional differences in bottom-up indicators and possible drivers to be integrated into spatial management policies (Tweddle et al. 2018).

Bottom-up indicators: top predators -gannet, harbour porpoise and harbour seal
Gannet, harbour porpoise and harbour seal appear to be the identified top predator indicators of bottom-up effects.Gannets are the most generalist foragers of the seabird species, foraging at a broad depth range within the water column (Mitchell et al. 2020), and harbour porpoise are the only marine mobile in this study not to exhibit centrally placed foraging, with known large spatial distribution shifts (Evans and Waggitt 2020) and harbour seals take a wide variety of prey which can vary seasonally and across regions (SCOS 2021).For these reasons we speculate that species with more flexible diet and/or locational needs are better bottom-up indicator species.Seabirds and mammals have long been recognised as being useful (more easily detectable) 'indicators' of the health of marine ecosystems (Parsons et al. 2008, SMRU 2016), with some specialists such as kittiwakes suggested as being the better indicators of environmental change (Wanless et al. 2007); however, results from our study question if we have been going about this the wrong way.It may be better to know what precisely the particular visible top predator species is good at indicating in terms of physical or 16000587, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ecog.06925 by Battelle Memorial Institute, Wiley Online Library on [19/12/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License biophysical ecosystem changes, as this work has showed for gannet, harbour porpoise and harbour seal to be indicators of bottom-up effects with numerous strong interactions with specific physical and/or biological ecosystem components.This increase in knowledge about what top predators are telling us about how the ecosystem is functioning and responding to environmental change is essential for effective spatial management, and we need to pay much closer attention to what indicator species are telling us about how specific components of the ecosystem are likely to be influenced by different climatic conditions or changed physical conditions due to ORE developments.

Top-down indicators: sandeel larvae and fish assemblages
As fishing has been demonstrated as a long-term dominant external stressor (Lynam et al. 2017), it is not surprising that all fish species/assemblages were identified as indicators of top-down fisheries driven forcing with changes in population trends that were dependent on the levels of fishing across the regions.For example, in the shallow central North Sea, a region characterised by stronger top-down fishing effects, the demersal assemblage was consistently predicted with a large increasing trend in response to scenarios with low fishing levels, i.e. any reduction in fishing pressure is positive for demersal fish.In contrast, in a less top-down, more bottom-up dominated region (Shetland/Orkney) neither low nor medium levels of fishing produced much of an increasing trend for either sandeel larvae, pelagic or demersal fish assemblages (Supporting information).For sandeels, much of this area has been closed for fishing for decades (Régnier et al. 2017), thus reducing fishing alone has little effect.These results highlight that we can use sandeel larvae and fish assemblages to assess ecosystem resilience to anthropogenic changes to identify what habitats (e.g.shallow central North Sea) and species (e.g.demersal assemblage) are more at risk.This approach will help to support effective spatial management of exploited species and avoid potential knock-on effects of ORE developments, such as fisheries displacement.

Bottom-up and top-down indicators: top predators -grey seal, guillemot and kittiwake
The results highlighted the following top predator species: grey seal, guillemot and kittiwake as indicators of assessing effects of not only physical effects of climate change and ORE developments, but also effects of fishing.Strong regional differences in colony demographics and population dynamics have been reported for grey seals which has made it challenging to identify the drivers of population change (SCOS 2021).Our results, based on the use of regional productivity data, have concluded that by identifying strong and consistent relationships, but with very different indicators in the different regions (Fig. 6a and grey seals plots in Supporting information), both physical (e.g.PEA) and biological (e.g.sandeel larvae) ecosystem components indeed drive production trends.Other modelling studies have discussed the role of grey seals in regulating the equilibrium of the ecosystem by exerting top-down control on their prey (Serpetti et al. 2017).
A very recent modelling study, based on assumed relationships rather than population trends, showed that the number of grey seals modelled tend to vary based on the interactive effects of both warming and fishing through changes in lower trophic levels and pelagic and demersal fish (Thorpe et al. 2022).With this range of new evidence, it is important that the spatial location of grey seals should be considered when it is being used as an indicator species.
Our findings have demonstrated that guillemot is also an indicator of predicting how the ecosystem is responding to the interactive effects of multiple bottom-up and top-down changes, which is in line with previous findings about seabirds, highlighting the complex ways in which climate interacts with other drivers with impacts on prey availability (Mitchell et al. 2020).Our results showed that demersal/pelagic fish assemblages (Fig. 6c and guillemot plots in Supporting information) drive guillemot breeding success in all regions, which possibly indicates the potential for guillemot to also be a strong indicator for top-down fishing effects.The time lags that we encountered in some of the annual predictions for guillemot (Fig. 6c) are, we suggest, due to structural uncertainties in either the species-specific and/or colony-specific drivers and whether their relationships with guillemot annual breeding success should be modelled inter-annually (i.e. between two or more years) versus intra-annually (e.g.within the same year), as the few studies on other species have shown inconsistent patterns between different colonies and their drivers (Carroll et al. 2015, Eerkes-Medrano et al. 2017).To address this, we will be investigating ecosystem models on a finer spatial scale with spatially and temporally rich baseline data (e.g.Isle of May), such that the understanding of mechanisms is achieved and can be subsequently transferred across larger regional and shelf-wide scales.
For kittiwake, the HDBN model did well in reproducing the observed trends in annual breeding success (e.g. decline in early 2000s and later increases, shown in Fig. 6e, Supporting information), suggesting that bottom-up forcing is mediated through interactions with primary production and zooplankton but also indirect influence from fishing is mediated through interactions with sandeel larvae.However, our results for kittiwakes failed to predict any explicit trends in response to the scenarios, which we explain is due to the complex effects across and between trophic levels, resulting from the interacting effects from other factors (e.g.fishing and temperature, and phenological effects, such as the timing of phytoplankton blooms) which might be more influential drivers governing this species dynamics (Carroll et al. 2015, Eerkes-Medrano et al. 2017).The ambiguity around the potential drivers for kittiwake is a good example reminding us that the strength of relationships with indicators that control population trends/breeding success can rise and fall over time (Trifonova et al. 2021); and, as this study shows, specifically for the deep central North Sea, so as the influence of bottom-up versus top-down effects can change over time.Therefore, more detailed understanding across a range of spatial (< 1 km through to 1000 km) and temporal (days through to years) scales is needed to pinpoint the mechanisms 16000587, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ecog.06925 by Battelle Memorial Institute, Wiley Online Library on [19/12/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License behind species dynamics, as studies have shown inconsistent patterns between different colonies and their drivers (Carroll et al. 2015, Eerkes-Medrano et al. 2017).Still, our results highlight kittiwakes as indicators of both bottom-up and top-down effects.Given the large declines in breeding abundance of this species (JNCC 2021), ongoing scrutiny is needed across habitats to continue to improve understanding of how both types of forcing can alter ecosystem dynamics via kittiwakes' relationships with other ecosystem components under future climate and anthropogenic change.

Conclusion
The effects of physical and biophysical interactions on populations and ecosystems, and how these vary with climate change and anthropogenic transformations, need to be better understood.Given the complexity of marine ecosystems, full understanding might be challenging, so we need to have pragmatic methods to make advances.Through the applied BN approach, we were able to make tractable predictions of the true dynamic nature of bottom-up versus top-down effects and their patterns across trophic levels and regional habitats, and their changes over time.This increases knowledge necessary to add to the traditional use of top predator population dynamics as separate aspects of marine systems, and will reduce uncertainties of the level of direct and indirect effects on populations across a range of trophic levels.Our results indicated that the ecosystem dynamics in some regions (Shetland/Orkney, west of Scotland) might be in a more desirable ('healthy') state (given the number of predicted increasing population trends).However, due to the interacting effects of climate (temperature, PEA and plankton production), fishing and future implications of ORE developments -including potential fisheries displacement -all these factors should be considered within current management and decision making.We have showed that the ecosystem dynamics in the shallow central North Sea might not be as robust to climate and/or anthropogenic changes; therefore, many trade-offs will need to be weighed up rapidly for the future sustainable management of marine ecosystems between different uses of our shallow shelf seas under the influence of future changes.
We also showed that the deep central North Sea region seems somewhat unique with both types of forcing (bottomup and top-down) leading to complex patterns of control on the ecosystem, suggesting intense scrutiny is needed within this region, given the future implications of hundreds of GW of ORE developments by 2050 (Wind Europe 2022).The identified indicators of bottom-up (e.g.gannet) and topdown effects (e.g.sandeel larvae), as well as indicators of both bottom-up and top-down forcing (e.g.grey seal) should be used to better guide strategic and integrated research approaches to monitoring and field surveys.Specifically, bottom-up indicators (e.g.NetPP and zooplankton) should be prioritized as they affect the availability of fulcrum pelagic fish species (Cury et al. 2000) that are the main prey of mobile top predators and larger commercial fish species.It has been shown that top predator foraging locations are due to fish availability being tied to locations of new primary production, including tidal fronts and patchy areas of subsurface chlorophyll -areas where fish are actively foraging in space (Scott et al. 2010, Embling et al. 2013, Scales et al. 2014).However, it is also necessary to include top-down indicators (e.g.sandeel larvae and fish assemblages) to be able to predict the changes in fish availability from the changes in fishing fleet distributions and levels of catch induced by the changes in fishing levels due to ORE developments.Finally, insights on the identified indicators of both bottom-up and top-down driven effects, with better understanding of what type of ecosystem change (physical or biophysical) they are indicating and their predicted responses to combined climate and anthropogenic change, should also be used to evaluate accurate cumulative effects through impact assessments to support MSP and to enable evaluation of trade-offs to provide the most sustainable future spatial use of our seas.

Figure 1 .
Figure 1.The spatial boundaries of the four regions: Shetland/Orkney, west of Scotland, deep and shallow central North Sea.Map based on Fig. 1 in Trifonova et al. (2021).

Figure 2 .
Figure 2. (a) General structural form of the dynamic Bayesian network model with a hidden variable (HDBN) where X 1 …X N represents the set of variables and arrows denote conditional independence relationships.(b) The strongest data-driven relationships, identified in Trifonova et al. (2021), that were used to build the HDBN model for the Shetland/Orkney region in this study.Blue-coloured links indicate relationships with the physical indicators, red with fisheries catch, green with primary production components, orange with zooplankton assemblages and purple with sandeel larvae and fish assemblages.Four different zooplankton assemblages were used to develop the models in this work; however, for visual purposes, relationships with only A4 are shown.Symbols used to denote the ecosystem components are shown below the relationships.

Figure 3 .
Figure 3. Baseline predictions (blue line) versus real data (black line) and moderate (green line) and extreme (red line) scenario predictions generated by the Bayesian network model with a hidden variable (HDBN) model for the bottom-up lower trophic level indicators and their changes over time.Solid black symbols represent the indicators used to build the relationships for the modelling network structures, whilst coloured symbols indicate scenario changes driving those indicators.

Figure 4 .
Figure 4. Baseline predictions (blue line) versus real data (black line) and moderate (green line) and extreme (red line) scenario predictions generated by the Bayesian network model with a hidden variable (HDBN) model for the bottom-up top predator indicators and their changes over time.Solid black symbols represent the indicators used to build the relationships for the modelling network structures, whilst coloured symbols indicate scenario changes driving those indicators.

Figure 5 .
Figure 5. Baseline predictions (blue line) versus real data (black line) and moderate (green line) and extreme (red line) scenario predictions generated by the Bayesian network model with a hidden variable (HDBN) model for the top-down indicators and their changes over time.Solid black symbols represent the indicators used to build the relationships for the modelling network structures, whilst coloured symbols indicate scenario changes driving those indicators.

Figure 6 .
Figure 6.Baseline predictions (blue line) versus real data (black line) and moderate (green line) and extreme (red line) scenario predictions generated by the Bayesian network model with a hidden variable (HDBN) model for both bottom-up and top-down indicators and their changes over time.Solid black symbols represent the indicators used to build the relationships for the modelling network structures, whilst coloured symbols indicate scenario changes driving those indicators.

Table 1 .
Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ecog.06925 by Battelle Memorial Institute, Wiley Online Library on [19/12/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Summary of data.See Supporting information for more detailed information on the data sources and citations.

Table 2 .
Sum of squared error (SSE) of ecosystem components predictions generated by the Bayesian network model with hidden variable (HDBN) models.Values highlighted in bold are the most accurately predicted (values of SSEs that are less than 25.00) ecosystem components per region.The symbols denote component-specific interactions that are used to build the HDBN models across the four regions.PEL, pelagic fish assemblage; DEM, demersal fish assemblage.