A Bayesian network to simulate macroinvertebrate responses to multiple stressors in lowland streams

Aquatic ecosystems are affected by multiple environmental stressors across spatial and temporal scales. Yet the nature of stressor interactions and stressor-response relationships is still poorly understood. This hampers the selection of appropriate restoration measures. Hence, there is a need to understand how ecosystems respond to multiple stressors and to unravel the combined effects of the individual stressors on the ecological status of waterbodies. Models may be used to relate responses of ecosystems to environmental changes as well as to restoration measures and thus provide valuable tools for water management. Therefore, we aimed to develop and test a Bayesian Network (BN) for simulating the responses of stream macroinvertebrates to multiple stressors. Although the predictive performance may be further improved, the developed model was shown to be suitable for scenario analyses. For the selected lowland streams, an increase in macroinvertebrate-based ecological quality (EQR) was predicted for scenarios where the streams were relieved from single and multiple stressors. Especially a combination of measures increasing flow velocity and enhancing the cover of coarse particulate organic matter showed a significant increase in EQR compared to current conditions. The use of BNs was shown to be a promising avenue for scenario analyses in stream restoration management. BNs have the capacity for clear visual communication of model dependencies and the uncertainty associated with input data and results and allow the combination of multiple types of knowledge about stressor-effect relations. Still, to make predictions more robust, a deeper understanding of stressor interactions is required to parametrize model relations. Also, sufficient training data should be available for the water type of interest. Yet, the application of BNs may now already help to unravel the contribution of individual stressors to the combined effect on the ecological quality of water bodies, which in turn may aid the selection of appropriate restoration measures that lead to the desired improvements in macroinvertebrate-based ecological quality. © 2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
The ecological status of water bodies is affected by multiple stressors acting over multiple spatial and temporal scales (Allan et al., 1997;Frissell et al., 1986;Roth et al., 1996), such as increasing water temperature, changes in flow, reduction of morphological heterogeneity and increasing nutrient loads ( Friberg, 2010 ;Tockner et al., 2010 ). Moreover, these stressors may interact, having synergistic, antagonistic or additive effects on the ecological sta-tus of freshwater bodies ( Jackson et al., 2016 ;Piggott et al., 2015). The combined effects of these multiple interacting stressors are, however, still poorly understood ( Folt et al., 1999 ;Jackson et al., 2016 ), since stressor-response relationships observed in controlled experiments are specific to organisms, stressors and environments and are therefore difficult to extrapolate to the field ( Jackson et al., 2016 ).
The lack of understanding of the combined effects of multiple interacting stressors may also explain why knowledge of the effect of specific management interventions on ecological water quality is still limited ( Palmer et al., 2005 ;Pander and Geist, 2013 ). Consequently, a high proportion of restoration measures are ineffective, even now ( dos Reis Oliveira et al., 2020 ;Palmer et al., 2010 ).
Hence, to increase the effectiveness of restoration measures, we first need to increase our knowledge of how ecosystems respond to multiple stressors and to unravel the contribution of the individual stressors to their combined effect on the ecological status of waterbodies.
Model simulations provide the opportunity to relate the state of an ecosystem to environmental changes as well as to restoration measures and simultaneously provide understanding of the underlying ecological interactions. Consequently, models may be used to predict the effects of management interventions on ecosystem states in space and time and thus provide valuable tools for water management. Over the last decades, several ecological prediction models have been developed, ranging from mechanistic representations of environmental processes to food web models and statistical data-driven models ( Janssen et al., 2015 ). For the latter, techniques have been used such as decision trees, artificial neural networks, generalised linear and additive models, fuzzy logic models and Bayesian Networks (BNs) ( Pistocchi, 2018 ). The construction of such statistical ecological prediction models can be data-driven, knowledge-based or a combination of both ( Mouton et al., 2009 ;van Echelpoel, 2020 ). A review of the advantages and drawbacks of selected modelling techniques indicated that BNs are promising tools for the combined application of expert knowledge and ecosystem measurements ( de Vries et al., 2020a ;van Echelpoel, 2020 ).
BNs are causal network models in which nodes depict (environmental) factors and in which dependencies between nodes are expressed as probabilistic relationships . The main advantage of this type of model is that the full range of available knowledge on cause-effect relations can be used, originating from experts, mechanistic modelling output, literature and experimental and observed data ( Landuyt et al., 2013 ;McCann et al., 2006 ), integrating the scattered knowledge on cause-effect relations in water bodies. Moreover, in these models, samples with incomplete datasets can still provide knowledge ( Barton et al., 2012 ). The uncertainty associated with the input data is explicitly accounted for, and the predicted outcome is reported as likelihoods ( Uusitalo, 2007 ). In addition, BNs provide a visualisation of the causal relationship between the predictors, which helps with communication of the model. Limitations of this model type are the lack of representation of feedback-loops, and the requirement for discretising continuous data ( Uusitalo, 2007 ). However, for evaluating stressor-effect relations in water bodies we considered that the numerous advantages of BNs outweighed these drawbacks. The aim of the present study was therefore to develop and test a Bayesian Network for simulating responses of stream macroinvertebrates to multiple stressors. Since we anticipated that stressor-effect relationships would be context-specific, model predictions were approached using water-type and region-specific relationships ( de Vries et al., 2020b ). To this end, a BN model was developed which included the links between macroinvertebratebased ecological quality and stream characteristics for a single water type, the temperate, sandy lowland streams within the Northwestern European plain. The availability of an extensive dataset with measurements of multiple stressors and ecological responses for Dutch lowland streams enabled us to develop this BN-model. The developed model was then applied to predict the influence of stream restoration management scenarios on ecological quality as represented by macroinvertebrates.

Study area
The studied lowland streams were located on the ice-pushed ridges in the Veluwe area in the centre of the Netherlands ( Fig. 1 ).
The land use in the catchments consisted of agricultural fields, urban areas and deciduous and coniferous woodlands. Mean annual rainfall in the study area was 850 mm and daily temperature varied between -16 and 29 °C. The flow velocities in the groundwaterand precipitation-fed streams varied strongly (10-80 cm/s). The stream bottoms consisted of sand and gravel.

Macroinvertebrate and environmental data
The data was collected by the Dutch Water Authority 'Vallei and Veluwe' during regular monitoring programmes over the period 1981-2017. In total, 208 sites in the upper courses of the lowland streams were selected. At these sites macroinvertebrate abundance data was collected as a part of regular monitoring programs. For each macroinvertebrate sample the ecological quality ratio (EQR) was calculated according to the Dutch assessment system, which expresses the ecological quality of a water body (ranging from 0-1.0) as a fraction of the reference situation (1.0) ( Van der Molen et al., 2016 ). In addition, for each macroinvertebrate sample, the mean preference score (ranging from 1-5) for several environmental variables of all species present in that sample was calculated using relative abundance frequencies. To this end, an environmental preference dataset was used ( Verberk et al., 2012 ). Environmental variables that were monitored at the same locations and at the same moment as the macroinvertebrate samples included water temperature, dissolved oxygen concentration, stream velocity, shading, total phosphorous concentration, biological oxygen demand, chlorophyll concentration, stream gradient, silt cover, macrophyte cover, coarse particular organic matter cover, and the presence of wood and gravel. However, not all variables were measured at all sites and on all occasions and therefore only environmental monitoring data was included when macroinvertebrate abundance data and at least a single environmental variable were monitored simultaneously. This resulted in a set of 933 samples ( Fig. 1 ).

BN theory
In short, BNs consist of causal network structures in which nodes, representing important system variables, are related to each other through arrows, representing dependencies ( Charniak, 1991 ). The state of a node is determined by the states of its parent nodes. This approach is described by Bayes' theorem, in which prior probabilities are updated given the likelihood of the data to generate a posterior probability distribution ( Ellison, 2004 ) . The type of relation between a node and its parent nodes, as well as the associated uncertainty, are recorded in a conditional probability table (CPT).
CPTs can be based on multiple types of data, including expert knowledge, process-based modelling output, literature-based values, experimental and observational data. Observational data is preferred, but when gaps are present in the dataset, other types of evidence may be used to quantify the relationships between the nodes. CPT relations that are initially based on expert knowledge can also be trained by using field observations. Hence, BNs have the advantage of being able to deal with incomplete datasets, and of providing ways to combine different sources of knowledge ( Uusitalo, 2007 ).
A BN captures relations between a set of variables, which may be uncertain, probabilistic, or imprecise. When the predictions are used in decision making, the explicit reporting of the associated uncertainty and the variability in the model results provides an advantage of this approach over deterministic methods that lack this reporting . Another advantage of BNs is that calculations can be made in the two directions of the arrows between the nodes: the values of child nodes can be calculated given the values of the parent nodes and vice versa. Conse- quently, BNs can be used to predict the outcome of scenarios given a set of causal variables, but also as diagnostic tools to deduce the probabilities of causes given the observed consequences, by using the dependencies in the model structure backwards ( Barton et al., 2012 ;McCann et al., 2006 ).

Model development
A five-step model development process was adopted  following several guidelines ( Aguilera et al., 2011 ;Chen and Pollino, 2012 ;Landuyt et al., 2013 ). 1) Model structure: An influence diagram was set up showing the causal relations between the environmental variables and the macroinvertebratebased EQR. 2) Model parametrisation: The CPTs picturing the relationships between the nodes in the model structure were defined using expert knowledge, and model-based and literature-based relationships. In this step also continuous variables were discretized. 3) Model training: The CPTs were trained using observations. 4) Model testing: Model performance was tested using independent observations. 5) Model application: A final model version for application was trained using all available data.
Model structure (step 1) The constructed BN represents a causal network of the environmental factors that influence macroinvertebrate assemblages. The chosen outcome variable was the EQR. The development of the early stages of the model structure was described in Skeffington et al., 2014 s. Based on literature and input from stream macroinvertebrate experts, the key environmental factors that influence the EQR in this specific water type and region were selected, including temperature, oxygen concentration, flow velocity, food quality and substrate variability ( Sandin and Johnson, 2004 ;Verberk et al., 2012 ;Verdonschot et al., 1998 ). These factors score the response to the environmental variables on a scale from 0-1, thus giving the user the opportunity to see which stressor is most limiting for a high EQR. The optimal values of these key factors were based on water-type and region-specific preferences of the reference macroinvertebrate assemblage ( Verdonschot et al., 20 0 0 ). Next, predictors of those environmental variables that span local to regional spatial scales were included, resulting in the model structure depicted in Fig. 2 . The selection of the included variables was made to cover all relevant processes from reach-to catchment scale, but at the same time aimed to limit model complexity. In addition to the environmental variables, five nodes were included that represent the average preference score of the observed macroinvertebrate assemblage for the variables flow velocity, CPOM, silt, gravel and wood. These scores are based on a preference database that lists species-specific preferences for environmental variables, based on experimental and distribution data ( Verberk et al., 2012 ), and range from 1 (low preference) to 5 (high preference). Inclusion of these nodes served as an additional source of information about the quality of a water body based on the preferences of the local macroinvertebrate assemblage.
Model parametrisation (step 2) Initially, the relationships between the nodes in the model were based on literature, mechanistic modelling outcomes and expert judgement ( Table A.1 ). When possible, a relationship was expressed as an equation. Discretisation was either based on equal intervals or on equal frequency to assure an even spread of data ( Chen and Pollino, 2012 ). The number of discretisation classes for each node was either 4, 5 or 6, balancing a minimum resolution to picture environmental processes with sufficient data availability per class.
Training and testing of the network (step 3 and 4) Step 3) and 4) were combined in a k-fold cross-validation context ( Marcot, 2012 ). The dataset was split into 3 parts, of which 2 parts were used for training the network and 1 part for testing the network performance. For a stable validation outcome among trained model variations, subsets of data for cross-validation can be made using stratified classes. Although this was not applied here, the subsets of data did have a similar distribution of EQR classes. Model training and testing was done for 3 consecutive runs, where in each run another part of the dataset was used for testing. The resulting metrics were then averaged for overall model performance. Cross-validation was performed on several model variations to test the influence of differences in structure and parametrisation on the prediction performance, i.e. dis- cretisation method, number of discretisation classes and the inclusion/exclusion of nodes representing macroinvertebrate preference data ( Table 1 ). Sensitivity analysis of the network was performed to identify the factors that had the strongest influence on the target node. It was thought more valuable to compare the model results as a continuous EQR-value to observed continuous values, as the use of EQR classes would not be informative enough in practice. Conventional metrics such as the number of correctly classified instances only use the classified output of the model, expressed as discrete values, and are therefore less suitable for testing the performance of the model in predicting the actual continuous EQR. Therefore, in this study, performance was tracked using the correlation between the observed and predicted EQR-scores ( Marcot, 2012 ), although the performance might be more strictly assessed than by using class-based metrics. The EQR scores used were expected values predicted by the target node, which is the average of the discretized classes weighted by the probability of occurrence.
To develop and test the BN model, the Netica BN software was used as a modelling shell (Norsys, 1998). This software provides a graphical user interface, can handle input of continuous data, provides ways to perform sensitivity analysis and can work in batch mode to more easily run the model for multiple sites.
Scenario analysis (step 5) Based on the performance analysis, the best performing model variation was selected. Next, this model was trained with all available data ( Marcot, 2006 ) ( Table A.3 ) and subsequently applied to predict the influence of stream restoration management scenarios on macroinvertebrate-based ecological quality. To this end, sites  were relieved from either one or multiple stressors ( Table 2 ), in which only combinations of scenarios that targeted at least two different key factors were considered. A comparison was made between the effect of removing single and multiple stressors per stream or stream stretch.

Results
The first steps in developing the model were constructing the network structure and then informing it with literature-and expert-based knowledge (step 1 and 2). In the third step, the knowledge-informed network was trained using actual monitoring data. Training the network with monitoring data resulted in adjusted probabilities, directly affecting 17 nodes, with a maximum change in probabilities of 60%. For instance, there was a relatively high number of high-gradient streams in the dataset. Consequently, the corresponding probabilities were adjusted such that for any new datapoint without observations for stream gradient, the model assumed a higher prior probability that it is a high gradient stream. This in turn altered the prediction of the status of the target node, the ecological quality expressed by the EQR. This adjustment is reflected by the larger bar in the top right node of Fig. 3 b compared to Fig. 3 a.
In step 4 of the model development, the network was tested using a part of the dataset applying a 3-fold cross-validation. The performance of the tested model variations ( Table 3 ) was expressed as correlations between the observed and the predicted EQR scores and showed scores up to 0.35, expressing a relatively poor predictive performance. The model variation performing best was obtained by incorporating 6 discretisation classes, using equal frequency discretisation, and including the trait preferences nodes.
The sensitivity analysis shows the influence of the environmental factors on macroinvertebrate-based EQR, the target node, in decreasing order ( Fig. 4 ), revealing that temperature, velocity and CPOM had the strongest influence, whereas the factors macrophytes, stream gradient and total phosphorous concentration had a very limited influence. This best performing model variation was subsequently applied to compare the model EQR predictions with the observed EQR per stream ( Fig. 5 ), showing that in half of the cases, the EQR was predicted well.
This concerned mainly streams with a relatively low EQR. In contrast, in the other half of the streams, having a relatively high observed ecological water quality, the predicted EQR was lower than the observed EQR. In these underpredicted cases, either the model judged the environmental variables too severely, or the observed EQR was overrepresenting the actual ecological quality, as a result of a too optimistic underlying assessment system. To gain more insight into these underpredictions, the two most deviating cases, the Zwaanspreng and Egelbeek, were considered in more detail. To this end, the average preference score of the macroinvertebrate assemblage was calculated for each sample for the factor flow velocity, one of the most influential environmental variables ( Fig. 4 ) and the only factor for which enough preference data was available ( Fig. 6 ). The mean flow velocity preference scores for these samples was relatively high (mean preference score: 3.6 out of 5), and also the flow velocity was generally high (mean 0.31 m/s). Likewise, the observed EQR was also high (mean 0.81). This points to the model having underpredicted the EQR of these sites, based on the factor flow velocity. However, flow velocity is only one of the factors determining the EQR and therefore, also other factors might have contributed to the underestimated EQR. However, as there is insufficient data available for the actual assemblage preference for the other environmental factors, we could not evaluate the contribution of these factors to the underpredictions of the EQR.
In the scenario analysis (step 5), the model was used to predict the effect of relieving the stream from single and multiple stressors on the target node, the EQR. For all streams combined, significant differences in EQR were observed when compared with the current conditions ( Fig. 7 a). Yet, for scenarios involving the relief of a single stressor, more negative than positive effects were observed. In contrast, when a combination of stressors was removed, the majority of scenarios showed positive effects on the EQR. In some of the streams, the effects of taking away the stressors could not be predicted (not shown), which might be due to inconsistencies between the observed and the scenario-based variables in the model, where nodes receive contradicting input. When the scenario effects were considered per individual stream, there was a high variation in the results ( Fig. A1 ). Nevertheless, for half of the streams clear management effects were still observed. To illustrate this, the Hierdense beek and Tongerense beek were considered in more detail, because these streams showed the clearest effects of stressor relief and had the largest dataset, respectively. Moreover, the specificity of the predictions was increased when the samples were grouped in stream stretches that represent a specific waterbody subtype within similar surrounding conditions, as can be seen for the upstream stretch of the stream Hierdense beek ( Fig. 5 ). For this stretch of the Hierdense beek and for the stream Tongerense beek several positive effects of management scenarios were observed ( Fig. 7 b, 7 c). Relieving the stream from most single stressors and stressor combinations increased the EQR. In contrast, scenarios involving an increase in velocity showed negative effects on the EQR, except when this measure was combined with increased CPOM cover, where a strong positive effect was seen. Especially the combined approaches that increased CPOM cover and flow velocity or CPOM cover and substrate quality (wood and gravel presence) had a positive effect on the mean EQR of both streams.

Discussion
The aim of our study was to develop and test a Bayesian Network for simulating macroinvertebrate-based ecological water quality based on the responses of stream macroinvertebrates to multiple stressors. The model was developed for a specific water type in a single region, where multiple stressors affected the stream ecosystem quality. Although surrounded by substantial margins of uncertainty (as seen in Fig. 7 ), the BN clearly showed the positive influence of restoration measures on the ecological quality of the studied lowland streams. Below we will discuss the performance of the BN in the scenario analyses, the complexity of 6 Fig. 4. Sensitivity analysis for the best performing model variation. For abbreviations, see Table 2 . predicting the effects of multiple stressors on macroinvertebratebased ecological quality and how this approach can be applied in stream restoration management.

BN model development
BNs are promising as tools in restoration management, as they offer a way to include expert knowledge with associated uncer- Fig. 5. Comparison of mean EQR model predictions and mean EQR observations for the studied streams. Boxes are inter-quartile ranges (IQR, 25 th to 75 th percentile) with whiskers extending to -/ + 1.5 * IQR. Statistical pairwise differences were calculated using Wilcoxon test, * : p < = 0.05, * * : p < = 0.01, * * * : p < = 0.001, * * * * : p < = 0.0 0 01.  Table 2 . Asterisks indicate significant differences in EQR compared to the current conditions (CC) (Wilcoxon test). a) Model predictions of restorations for all streams combined (based on a total number of 9891 model runs, with a variable number of model runs per scenario due to inconsistencies). b) Model predictions of restoration scenarios for the Hierdense beek (upstream) (based on a total of 215 model runs). c) Model predictions of restoration scenarios for the Tongerense beek (based on a total number of 1552 model runs). tainties in combination with monitoring data and output from process-based modelling. In addition, the explicit expression of the associated uncertainty and the clear visualisation of the causal model structure are advantages of this technique, which make it a suitable tool to support water managers in decision making ( Barton et al., 2012 ;Uusitalo, 2007 ). However, to be able to train the knowledge-informed network using observations, for each pair of connected nodes, there should be a set of data available that covers the full range of all possible combinations of node states ( Cain, 2001 ). This is a requirement which may be difficult to comply with in practice, especially when the focus is on a specific water type.
In the development of our BN model, choices had to be made to deal with the inherent complexity of aquatic ecosystems. To this end, the main predictors of the macroinvertebratebased EQR were selected. However, there was a trade-off between the desired model complexity and the availability of training data, where a lack of data would decrease model performance (see also Marcot et al., 2006 ). Hence, to select the optimal model variation, multiple model structures and parameterisations were tested. In the comparison of these slightly adjusted models, equal frequency discretisation gave better predictive performance than discretisation based on predefined class boundaries (however, compare Boets et al., 2015 ). Although all discretisation methods imply a simplification of continuous data ( Aguilera et al., 2011 ), using equal frequency discretisation ensures that each class of a node is represented equally in the data, which supports a better training of the conditional probabilities in the network. Also, the inclusion of the preference nodes slightly improved the predictions.
Our results showed the impact of restoration measures in the scenario analysis compared to the current situation, but the overall absolute performance of the BN model was still limited. Especially for streams with a high observed ecological quality, the EQR was underpredicted by the model. This might be partly due to gaps in the dataset, consequently, it was not possible to train each knowledge-based CPT with observed data ( Table A.3 ). A possible explanation may also be that the model only predicts the effects of changes in environmental factors on macroinvertebrates, whereas in reality, dispersal and biotic interaction filters also determine the macroinvertebrate assemblage composition and therefore the ecological quality of a specific site ( Poff, 1997 ). In addition, the interactions between environmental factors are not completely understood and cannot be fully incorporated. Moreover, the model is static and therefore assumes that the assemblage is in equilibrium with the environmental conditions at each site ( Austin, 2002 ), but this is not always the case ( Belyea and Lancaster, 1999 ;Wiens, 1984 ). Together, these complexities, which are not well understood, could have influenced the performance of the model.
Therefore, the current model is not yet thought to be acceptable for application as such, given that predictions are not yet in accordance with the observations. With a more complete dataset, testing of additional model variations and an increased insight in stressor interactions to improve model relations, this model type might be further applied as a tool in restoration management.

Multiple-stressor effects on macroinvertebrate-based ecological quality
Waterbodies are generally subjected to multiple stressors ( Birk, 2018 ). This creates a complex task for water managers who aim to improve the ecological status of stream ecosystems. In addition, stressor interactions may take place that either enhance the added effects of additional stressors (synergism) or decrease these combined effects (antagonism) ( Folt et al., 1999 ). Such interactive effects are specific to stressors, organisms and environments and consequently are difficult to predict ( Jackson et al., 2016 ). With this added complexity, simulating ecological quality remains a complex task, as our study showed.
Despite this complexity, the present BN scenario analyses showed that ecological quality can be improved when the streams are relieved from specific stressors or combinations thereof, whereas other restoration scenarios may prove to be less effective. For the studied streams, the strongest positive effect resulted from increasing flow velocity in combination with the presence of CPOM, whereas only improving flow velocity yielded no positive effect. This may be explained by the interaction between flow velocity and CPOM: only increasing stream velocity would not safeguard the variation in flow required for patches of coarse material to persist, providing necessary habitat for stream organisms ( de Brouwer et al., 2019 ). To gain more insight in the interaction between flow velocity and CPOM cover, these key environmental factors could be included in the network in more detail than they are now. The added value of combining restoration scenarios was also observed for the scenarios that enhanced the presence of wood and gravel substrate and the cover of CPOM. In the scenarios where the streams are relieved from the individual stressors, already a positive impact is seen, but when the stream is relieved from both stressors simultaneously, the macroinvertebrate-based ecological quality improves even more than expected based on the contributions of the single stressors, which could point at a positive synergistic interaction of stressor relief. Similarly, such interactive effects on macroinvertebrates have been observed for other environmental stressors (Beermann et al., 2018;Jackson et al., 2016 ). Yet, for the other combined scenarios in the present study where multiple stressors were adjusted, no interactive effects were observed. Indeed, also for other water bodies it was reported that additive effects of multiple stressors prevail (Gieswein et al., 2017). However, the current model is partly knowledge-informed and not completely based on data. This is especially the case for the target node, where the relationship picturing the combination of multiple stressors into a combined response was based on expert knowledge.
To better quantify the interactions between the stressors of interest, additional statistical analyses could be carried out on a more extensive dataset ( Feld et al., 2016 ;Glendell et al., 2019 ). In addition, experiments may help to disentangle the interactive effects of multiple stressors ( Elbrecht et al., 2016 ;Verberk et al., 2016 ). Only when we have more knowledge about the nature and interactions of stressor-response relationships for specific species and complete assemblages, can we develop modelling of multiple stressor impacts further. In turn, the application of models can show us where these knowledge gaps persist and where additional experiments might be needed to better understand underlying processes. Consequently, BNs and other approaches are complementary in their contribution to an increase of the understanding of multiple stressor effects.
Apart from the interaction between stressors, other studies showed that choosing measures based on identifying multiple stressors covering the entire catchment proved to be more effective ( dos Reis Oliveira et al., 2020 ;Feld et al., 2011 ;van Puijenbroek et al., 2019 ). This illustrates the significance of simultaneously considering multiple stressors over multiple scales for effective stream restoration.
In conclusion, the present model exercise demonstrated that applying different scenarios enhances the understanding of the effects of combinations of measures on macroinvertebrate-based ecological quality and may aid in selecting and prioritizing the most promising restoration measures, as discussed below. 9

BNs as tools in restoration management
Nowadays, in the practice of stream restoration, the selected measures are still strongly based on assumptions rather than proofs of their positive effects on the ecological status of stream ecosystem ( dos Reis Oliveira et al., 2020 ), which was the main motivation to perform the present study. The application of the current BN model indeed enhanced the insights into the possible effects of management scenarios on the ecological quality as represented by macroinvertebrates. Hence, these model predictions may be used to inform water managers which measures to prioritize in restoration management to effectively alleviate stress. This increases the chance that the applied restoration measures do indeed lead to the desired improvement in the ecological status of stream ecosystems. Whereas we used the BN model to show the relative impact of (combinations of) restoration measures on macroinvertebrate-based ecological quality, such models can also be used 'backwards' in a diagnostic approach to find causes for observed symptoms ( Feld et al., 2020 ;Trigg et al., 20 0 0 ). Ideally, a model performing well, tailored to the study area, would give insight into which environmental factors would produce most effect. For the manager, the next step would be to identify how these variables might be targeted, by linking these to actual restoration measures. However, this prioritisation is often not just based on the outcome of the model. In these scenario analyses, the use of site-specific knowledge would permit the manager to decide which variables to prioritize, for example, knowledge of the possibility and cost of certain measures, and of restoration efforts and disturbances that have taken place in the past.
As shown here, scenario analyses can be especially informative in situations where multiple stressors are acting. In the case of a single dominant stressor, a specific measure may be more easily selected, but for a situation with a more even contribution of multiple stressors, selecting and prioritizing restoration measures may not be straightforward. In these cases, scenario analyses may help to choose a combination of measures to alleviate the pressure on the ecosystem and to improve the ecological water quality.
The current model was designed and trained for a single area and water type. When applied to other areas, the main model structure can still be used as a starting point, although the choice for key environmental factors, parametrisation of the CPTs, calibration and validation should be carried out in a way tailored to the water type and region of interest.
Ultimately, in the application of BNs, challenges remain with the abovementioned complexities. In addition, Kaikkonen et al. (2021) list the remaining challenges of BNs used in environmental management, such as models lacking validation, unclear discretisation methods, and lack of clarity about the source of expert knowledge. Indeed, most BN applications fail to test the predictive ability of the model (Death et al., 2015). As described here, discretisation and validation of the model outcomes is not straightforward. Better reporting of such challenges associated with these technical aspects may therefore improve future robustness of BN applications. In addition, recent technical developments might further increase the possibilities of BN applications, such as the use of hybrid networks that can represent continuous variables without the information loss associated with discretisation ( Kaikkonen et al., 2021 ).
The success of applying BNs for similar purposes in the future depends on the availability of high-quality data and the possibility to include a more fundamental understanding of the complexity of ecosystems, with context-specific knowledge on how interactions between multiple stressors affect macroinvertebrate assemblages. The current approach has contributed to an increased understanding of the complexity of these aquatic ecosystems. Moreover, our study showed how BNs can be used in a scenario analysis to select and prioritize the most promising restoration measures.

Conclusions
In this study, the application of BNs for simulating the effects of multiple stressors on macroinvertebrate-based ecological water quality was tested. Although the predictive performance can be further improved, our application illustrated how these models can be used to increase our knowledge of how ecosystems respond to multiple stressors. To make predictions more robust, a deeper understanding of stressor interactions is required. Also, sufficient training data should be available for the water type of interest. Still, BNs allow us to make steps in unravelling the contribution of the individual stressors to their combined effect on the ecological quality of water bodies. This in turn may aid the selection of appropriate restoration measures that lead to the desired improvements in ecological water quality.

Author contributions
RAS, AJW and PFMV developed the initial BN. JdV further developed the BN and analysed the data. JdV, PFMV and MHSK led the writing of the manuscript. All authors contributed to the review and revision of the model and the manuscript, and approved the final article.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.