From responses of macroinvertebrate metrics to the definition of reference metrics and stressor threshold values

The main obstacles to monitoring aquatic ecosystems have always been the lack of data, the complex socio-economic context and the lack of specialised expertise that characterise the West African region. To overcome these challenges, it is necessary to invest in the definition of context-appropriate approaches that allow good ecological assessment, improve understanding of the functioning of these ecosystems and facilitate data collection. This study focused on using benthic macrofauna to assess the risks of moving away from Good Ecological Status towards the functioning of an anthropised system (Nokoué-Benin), based on the definition of reference values for macroinvertebrate metrics, stress thresholds and the responses of selected metrics to stressors. The approach used is a combination of a joint species distribution model and Bayesian networks. We used JSDM to select relevant metrics and generate posterior probabilities. We then converted these posterior probabilities into posterior response probabilities for each of the stress levels and fed them into a Bayesian network. We used the Bayesian network response predictions to define the reference values of the metrics and the stress thresholds derived from the probability density plots for the low pressure levels. An application of this approach was then carried out on a lagoon sampled during high and low water periods for three years, with 33 macroinvertebrate taxa present in all seasons and sampling points, and measurements of 14 environmental parameters used as application data. The relevance of the results, despite the small sample size, supports wider application of the approach in West Africa.


Introduction
Coastal lagoons are known for their great biological diversity and for the many ecosystem services they provide (WRI 2001).They are also ecotones with a fragile balance (Lloret et al. 2008) highly impacted by multiple human activities (Parry et al. 2007).Some of the anthropogenic impacts on the lagoons range from intense pollution to eutrophication and sedimentation, to the scarcity of water exchanges with the sea (Carlier et al. 2008;Souchu et al. 2010).In addition the threats extend to the effects of climate change, and often to the expansion of non-native and sometimes invasive species (Chapman 2012;Reizopoulou et al. 1996;Suarez-Menendez et al. 2020).
In Africa, the extent of impacts on lagoons is often unknown given the lack of monitoring data (Ahoutou et al. 2021).However, as population growth is likely to continue to increase strongly (UNECA 2016), the impact of multiple anthropogenic pressures is expected to be even more important in these tropical environments than in temperate regions (Ahoutou et al. 2021), while the use of these ecosystems is of paramount importance, as demonstrated in North Africa (El Mahrad et al. 2020).
Considering the distribution of organisms to account for the impacts of environmental stressors on the ecosystems and their ecological integrity (Carignan and Villard 2002;Siddig et al. 2016) is a very common practice (Birk et al. 2012).Among the biological components used as biological indicators, macroinvertebrates are often considered for their ability to respond to abrupt changes in quality and their sensitivity to pollution (Tutu 2017;Uherek and Pinto Gouveia 2014;Villantes 2015).
This paper presents a methodological approach aiming to assess the risks to ecosystem functioning in anthropized environments, by defining reference values of macroinvertebrate metrics as well as stress thresholds beyond which its trajectory would move away from good ecological status.We used Joint Species Distribution Models (JSDM) to simultaneously consider environmental factors and biotic interactions (Tikhonov et al. 2017) in the responses of metrics to stressors.Bayesian Networks (BN), are used for their effectiveness in modeling biotic interactions (Staniczenko et al. 2017) and simulating the responses of macroinvertebrates to anthropogenic pressures (de Vries et al. 2021).BNs make it easier to consider synergistic, antagonistic, and additive interactions in the pressures/impact's relationships.Indeed, incorporating the effect of interactions in the assessment of ecological status and in the establishment of so-called ''reference'' values is perceived as necessary (Miguet et al. 2019;Piggott et al. 2015), to the extent that its abstraction can lead to surprising responses to disturbances (Paine et al. 1998).Our approach included 3 steps (i): Identifying the relevant stable metrics that best respond to the stressors (physicochemical parameters, artificial substrates, and non-natural land cover of the watershed); (ii) Predicting the responses of the metrics to multiple stressors simultaneously and then (iii) Deducing the reference and stress thresholds values.The Nokoue ´, a lagoon in the Gulf of Guinea, subject to strong environmental constraints linked to various activities on the lake (fishing, acadjas…) and its watershed served as a framework for the development of the approach.According to previous monitoring work, its quality seems to be decreasing (Gnohossou 2006): a loss of 42% of its fish species richness, for example, was observed between the 1960s and 2000 (Djihouessi et al. 2019) and there is an urgent need for effective impact measurement.

Study area
The Nokoue ´(Fig.1), located between parallels 2°24 ' and 2°37' North and meridians 6°23' and 6°28' East (Mama et al. 2012), is a lagoon in the Gulf of Guinea in the Republic of Benin, West Africa.The Oue ´me ´River and the So ˆRiver in the North are its freshwater tributaries, while the Porto-Novo lagoon in the East, the Cotonou Channel and the Totche `Canal are its southern limits (Mama et al. 2012).

Macroinvertebrates data
The macroinvertebrate data were extracted from the PhD thesis of Gnohossou (2006).Two sampling campaigns were carried out during low water seasons (February 2005/February 2006) and two others during floods (September 2004/September 2006), on 80 stations (Fig. 1).The sampling consisted of collecting species present on a structure made of basket (diameter: 20 cm; height: 14 cm; mesh size: 0.5 cm), filled with bottom substrate and branches/fagots placed on the bottom for a one month immersion.This sample was completed by the use of a modified Surber of 50 cm width and 30 cm height with a 300 mm mesh net, dragged on the bottom in the manner of a dredge over a distance of about 3 m using a rope (Gnohossou 2006).
We described the 33 species based on seven traits: Life cycle duration, locomotion and substrate relation, food, feeding habits, substrate preference, trophic level and saprobity.The trait of each taxa was determined according to the literature (Bis and Usseglio-Polatera (2005), Demars et al. (2012) and Tachet et al. (2010) and completed by dedicated scientific and online resources (inpn.mnhn.fr/;www.marinespecies.org/;species-identification.org/; www.biotaxa.org/Zootaxa/).Some traits of taxa initially defined in Bis and Usseglio-Polatera (2005) were adapted according to information provided by experts of the lagoon.For some taxa, we used the experience of native fishermen as well as information from commercial journals.Each macroinvertebrate trait was subdivided into one or more secondary traits, and for each of these, three modalities were defined using fuzzy coding (Chevene et al. 1994;Tachet et al. 2010).In this coding system, for a given trait, a value of zero indicates the absence of affinity of a taxon and a value of five correspond to very high affinity.See Appendix SI-Tables [Tables 1 and 2] for lists of taxa and traits.For each sample point, diversity was calculated as the total number of taxa or traits observed.

Calculation of macro-invertebrate metrics and expected responses to pressure
Fifty-four candidate metrics related to the diversity, composition and abundance of species considered from a taxonomic or functional point of view (sharing a trait) were selected (SI-Appendix-Tables [Table 3]), for their effectiveness in accounting for environmental alterations in previous studies and their availability.Abundance metrics are the ratio of the total number of individuals constituting the taxa of a target group at a given sampling point to the total number of individuals of all taxa present at that sampling point.Functional metrics were obtained as the ratio of the total number of individuals sharing a trait for a given sample point divided by the total number of individuals for that sample point.A first group of candidate metrics included four species diversity indices: Simpson's k index (Sim_iL), Margalef's index (dMa), Berger Parker index (Bpi) and Shannon-Wiener diversity index (He).For these 4 diversity metrics, we expected lower overall species richness and diversity in the face of anthropic pressure (Erasmus et al. 2021).
A second group included metrics related to the taxonomic composition of the sample: the relative abundance of amphipods (Aamp), bivalves (Abiv), chironomids (Achi), oligochaetes and naididae (AOl_Nai), polychaetes (Apoly), mayflies (AEph) insects (AInsect), mollusks (Amol), gastropods, oligochaetes, and dipterans (AGOLD), mollusks and chironomids (Moch), gastropods (Agast), crustaceans (Acru), non-insect taxa (Ano_Insect), and potamids (Apot).For crustaceans, we would expect a decrease in abundance with sewage discharge, desalination (Navarro-Barranco et al. 2020a, b) and siltation, although they are poor indicators of low environmental disturbance.Bivalves are expected to increase in abundance due to their resistance to environmental stress (Helmholz et al. 2016); chironomids are associated with high disturbance (Lencioni et al. 2012); the same trend applies for annelids; oligochetes increase under organic pollution (Kazanci and Girgin 1998), while ephemeropteran abundance decreases Fig. 1 Reproduced map of the study area and sampling points; on the map, the ellipses indicate seven zones in which stations are distinguished.The station codes are formed by a letter designating the zone (E for mouth, V for Veˆki, G for Ganvie´, Z for Zogbo, C for Chenal and A for Atchonvi), then a letter designating the type of habitat (N for bare lake or open water, V for vegetation or macrophyte, A for acadjas) and, finally, the station number) (Gnohossou 2006) significantly under oxygen depletion (Menetrey et al. 2008); Coleopterans are associated with low phosphorus levels and are indicative of trophic status with negative responses to high anthropogenic stresses (S ˇidagyt_ e et al. 2013).

Pressure data and other environmental variables
The physico-chemical stressors considered in the analyses were BOD, pH, dissolved oxygen (DO), total phosphorus (TP), nitrates, ortho phosphates (Ortho_P), total kjeldahl nitrogen (TKN), dry organic matter (DOM), temperature (Temp), transparency (Trans) and salinity (Sal).These physico-chemical parameters were measured at the same time as the sampling of benthic macroinvertebrates (data from Gnohossou et al. (2006)).The hydromorphological alterations considered include: siltation relative to the alteration of the substrate, the depth (Depth) relative to the anthropic modification of the hydraulic geometry taken by sounding, the presence of acadjas (Acadjas) and absence of bank macrophytes (Macrophytes).Then an integrating stressor parameter represented by non-natural land cover (Oc_nN), calculated as the proportion of urban/agricultural cover in a 5 km circular buffer around a sampling point.Hydromorphological parameters and Oc_nN describe habitat degradation.The other forcing factors included as covariates were average wind speed (AWS), monthly average water level (MAWL), and year.All data on environmental and forcing factors came from published data (Hounye `me `and Daouda 2021).In addition, we included a station fixed effect (variable ''Sites'') and a random effect of the sampling units (Random sample).

Modeling
We used a Bayesian framework based on JSDM named HMSC (Ovaskainen et al. 2017).For all HMSC models run, the MCMC convergence needed to validate the reliability of the inference as well as the results is evaluated quantitatively using the criteria of effective sample sizes (ess) and potential scaling factors (psrf).Such that ess is very close to the theoretical value of the effective number of samples, while psrf is very close to one.The methodological approach consisted of four steps.

In a first step
To avoid the risk of multicollinearity between metrics and to save computational time, a redundancy analysis with linearity forcing between macroinvertebrate metrics was performed (step 1) using the ''redun'' function in the Hmisc-Harrell Miscellaneous package (Frank 2022).

Selection of relevant metrics-In a second step
The selected metrics were injected into a Bayesian framework (Drouineau et al. 2010).HMSC models enable Bayesian inference.The HMSC model was considered here as a pressure-impact model to simultaneously select the most robust and relevant measures that respond to anthropogenic alterations at the sampling site.A metric was selected if its model contained at least one stress parameter that was significantly expressed with or without natural covariates, such as the plotBeta graph (representing the b parameters, i.e. the responses of metrics to stress parameters), with at least 95% posterior probability of being positive or negative in the model indicating posterior support of the tracing threshold (supportLevel) greater than or equal to 0.7.In addition, the Pearson correlation coefficient (R 2 ) between expected and observed values had to be greater than or equal to 0.3 (Argillier et al. 2012) at the end of 3 successive HMSC analyses conducted on 2 training sub data sets corresponding to high and low flow periods, HMSC 1 and 2 respectively, and on a selection synthesis analysis sub data set performed on 1/3 of the stratified observations randomly taken on the whole data set (HMSC 3).With regard to the selection of stable metrics, the initiated HMSC 3 analysis allows more flexibility in the sampling period.Four chains (nChains = 4), 200 samples per chain (samples = 200), a long transient period (transient = 30,000), for a total of 90,000 iterations are used to fit the HMSC model with Bayesian inference.For all of the models, we specified a ''normal'' distribution.

The third step consisted of calculating response probabilities
The a posteriori probability of metrics responding to the different levels of a pressure was then collected from two equivalent HMSC models: HMSC_bru and HMSC_dis corresponding to models built and run respectively on the raw observation dataset and on the observation dataset manually discretized into different quality classes (Pressure Levels) for each of the pressure parameters.We assume that the a posteriori probability when a metric responds positively to a pressure (obtained via the output of the HMSC_bru model) is equal to the sum of the a posteriori probability that the same metric that has different levels of the same pressure, also responds positively (obtained from the output of the HMSC_dis model).It should be noted that, as in regression models with qualitative variables, the a posteriori probabilities obtained from the output of the HMSC_dis model all refer to a reference level (Pearson 2020).Taking the example of a pressure with 3 levels (High-Medium-Low) and whose reference modality is ''High'', the Prob1, Prob2, Prob3 and Prob4 values of the probabilities output from the HMSC_dis model can be represented as equations a (with Prob1 ?Prob2 = 1) and equations b (with Prob3 ?Prob4 = 1): The a posteriori probability of responses for each of the metrics knowing these different levels was determined from Eqs. 1 and 2 below, and from equations a and b, then normalized.The value of the intercession points between the negative and positive density curves, for a low-pressure level was considered as the reference threshold (Johnson et al. 2013).For this purpose, the a posteriori probability calculated in step 3 and the a priori probabilities directly derived from the interval discretized stressors (Table 1) applied to the data to transform them into discretized data, constituted the conditional probability table (CPT) describing the relationship between the nodes of the directed acyclic graphs (DAG) structure of the constructed BN.The DAG covers all relevant processes, with the aim of keeping the complexity of the model as low as possible.Indicators such as R 2 , Standard Deviation (SD) and RMSE allowed evaluation of DAG performance.The a posteriori probability of the calculated responses was compared to the expected responses of the metrics based on the literature to determine their validity.For unknown relationships, we kept those suggested by unsupervised learning.A better compromise, which would result from the different reference thresholds derived by this approach, would be to take the most severe reference value among all those obtained for the same stressor as the most appropriate for the definition of a reference station (discriminating stress thresholds).For analysis and interpretation, only the probabilities of metric responses to increased pressure are presented and discussed, due to the large number of pressure levels that make the number of derived posterior probabilities so large.Figure 2 shows the main steps of the implemented approach.

MCMC convergence diagnostics
Very low autocorrelation between consecutive samples were observed (Fig. 3).Similarly, effective sample sizes were very close to the number of effective samples, and potential scale reduction factors close to 1 for all HMSC models run (Fig. 3).This clearly indicates that the resulting posterior samples are representative of the true posterior distribution and inference reliability.

Selection of metrics (Step 1 1
Step 2) The HMSC models (HMSC 1, HMSC 2 & HMSC 3) selected seven metrics (Table 2) out of the 45 candidate metrics tested.This selection is the result of eliminating 26 metrics resulting from linearization and redundancy analysis, and finally 12 unstable metrics (R 2 \ 0.3).
During the flood period, the model (HMSC 1) selected 10 metrics, while 13 were selected from samples collected during the low flow period (HMSC 2).Among these metrics, eight were selected by both models.The best model predictions were obtained for Bpi, dMa, Sim_iL and Tpier.No seasonal differences were observed in the prediction quality of Bpi, Tperma and Tshred.The prediction quality of dMa, Sim_iL and especially Tpier is much better during high water than during low water.Conversely, Abiv and Tmud have a higher prediction quality during low water.Finally, HMSC 3 shows the instability of Tpier (R2 = 0.17).All selected metrics were sensitive to at least one stressor.
Three of the seven selected metrics are species diversity indices: Bpi (Berger-Parker index); dMa (Margalef index); Sim_iL (Simpson index k).Three others are related to species traits: Tmud (% substrate taxa with a preference for mud: Tanaidae, Peneaeidae, Portunidae and Caenidae); Tperma (% taxa permanently attached to the substrate: Balanidae, Mytilidae, Ostreidae and Serpulidae); Tshred (% shredding taxa: Aoridae, Melitidae, Amphilocus sp., Photidae, Munnidae, Dytiscidae, Noteridae) and one on the taxonomic composition of the sample: Abiv (abundance of bivalves (Corbulidae, Ostreidae)).None of the metrics selected in Table 2 showed any variation from year to year.While dMa was most sensitive to wind (AWS) and tides (MAWL), Bpi, Tmud and Tshred did not respond to any of the forcing factors.Tperma, Sim_iL and Abiv responded negatively to the effect of tides.Only Tpier eliminated by HMSC 3, responded positively to the effect of wind.

Metrics-pressure relationships
Trends in metric response probabilities (values [ 0.5) for high degradation levels shows many common responses between metrics but also some dissimilarities (Table 3).
The set of response probabilities for all other pressure levels is in SI-Appendix-Tables [Table 4].
As expected, Abiv, Bpi, dMa, Sim_iL, Tshred react negatively to the decrease in DO due to biodegradable organic matter, expressed by the parameter BOD, for which Tmud and Tperma react positively.The same is true for all parameters except Tmud, which reacted positively to an increase in the organic load (DOM).The response of Tmud seems a priori to be incoherent.The negative responses of Abiv, Bpi, dMa, Tshred and Tperma to hypoxia (DO) and to the stress factor TKN are in line with our expectations.Conversely, we observe a positive response of Sim_iL and Tmud to hypoxia.This is not what we expected.TP and orthophosphates are expected to increase, favouring Abiv and Tperma.This is not the case for Bpi and Sim_iL.The positive effect of Abiv and Tshred following the increase in nitrate and the negative response of Tmud, Tperma and Sim_iL seem to be as expected.This is not the case for Bpi and dMa.Acidification, however, induces the expected positive response of Sim_iL and  Tshred.As expected, increasing temperature reduces the taxonomic community of bivalves (Abiv), the most dominant taxa (Bpi) and those permanently attached to the substrate (Tperma).The expected responses were also observed at high salinities, which were unfavourable for the bivalve molluscs (Abiv) and Sim_iL and favourable for the other species.This is also the case for Abiv, Sim_iL and Tperma, which show a negative sensitivity to the decrease in transparency compared to the other metrics.
Concerning the hydromorphological stressors, siltation and acadjas positively influence Bpi, dMa and Tshred as expected.While acadjas have a positive effect on Tmud and Tperma, land artificialisation (Oc_nN) has a negative effect on the same metrics.This opposite response is also observed for these two metrics with respect to stressors such as depth and macrophytes.
The specification of the relationships between macroinvertebrate parameters and proxies, presented in the form of DAGs, clearly showed the influence of stressors on the metrics and thus on the ecological quality of the lagoon (Fig. 4).The DAG also showed coherent interactions between the different parameters (physico-chemical stressors, hydromorphological and forcing factors).Conversely, among the relationships in the model in Fig. 4, the one between 'Acadjas' and 'Macrophytes' does not seem obvious.
Regardless of the metric considered, the different values (Table 4) follow the same order of magnitude.The standard deviation between the thresholds obtained for the same stressor is small, except for DOM (3.55).With regard to the thresholds for severe stress, each of the different metrics, with the exception of Tperma, seems to be suitable for monitoring at least one stressor, since it has at least one discriminating threshold.This may indicate a complementarity between the metrics in the context of a water quality assessment.For salinity and according to the typology of transitional environments such as freshwater lagoons, we considered two thresholds.Thus, in the freshwater part, i.e. upstream, the threshold could be 2.08% and downstream, i.e. close to marine waters, would be closer to 8.86 %.

Selected metrics and coherence of metric responses
This study represents a significant advance in the ecological assessment of West African ecosystems, as all the selected metrics cover all facets of diversity except genetics, and we were able to determine responses to most of the identified stresses.Of these relevant metrics, only those, or combinations of metrics, that can distinguish between one or more quality classes or that have a strong correlation with the disturbance gradient should be included in the assessment (Hering et al. 2004;Schoolmaster et al. 2013).The approach implemented allowed the selection of 7 non-redundant, non-collinear, spatially stable metrics that could be calculated regardless of the season sampled and that responded simultaneously to different stressors.Three functional trait metrics 'Tperma', 'Tmud' and 'Tshred', known to respond well to anthropogenic disturbance in other environments (Vandewalle et al. 2010), were shown to respond to several stressors occurring in the Nokoue ´lagoon.The use of a %clingers metric similar to Tperma has already been described in the Overseas Assessment System to measure human impacts that alter food sources or habitat characteristics (Wagenhoff 2016).In addition, the bioindication potential of the taxa constituting ''Tperma'' has been justified by several authors: Balanidae (Poirrier and Partridge 1979), Mytilidae (Azizi et al. 2018), Ostreidae (Mai 2013), Serpulidae [Polychaetes] (Lu et al. 2017).
Tmud including crustaceans (Tanaidae-Peneaeidae-Portunidae) and insects (Caenidae) was included in a quality index assessing the trophic status of a New Zealand shallow estuary (Robertson et al. 2016).The inclusion of Tmud in this index allowed the identification of areas with high potential for the rapid accumulation of muddy sediments generally associated with increased organic matter and nutrient content and decreased oxygenation.The potential for indicator taxa included in this metric were demonstrated in the literature: Crustaceans (Tampo et al. 2021), Caenids [mayflies] (Dalu and Chauke 2019;Alhejoj et al. 2014).
Similarly, 'Tshred' which is only the association of crustaceans (Aoridae -Melitidae-Amphilocus sp.-Pothidae-Munnidae) and beetles (Dytiscidae-Noteridae) in our case, was selected as a good indicator of most of the stressors considered.This metric has previously been used to indicate changes in the processing of stream organic matter (Wagenhoff 2016).Crustaceans and Coleoptera that constitute this metric were demonstrated to respond well to several physico-chemical stress factors including pH, temperature, dissolved oxygen, salinity… (Kabore ´et al. 2016;Kang and King 2012;Correa-Araneda et al. 2010;Gesteira and Dauvin 2000;Navarro-Barranco et al. 2020a, b).However, despite the good indicator qualities of the 'Tshred' metric, intraspecific variation in among life history stages, seasonal type, or other ecological characteristics have been reported (Applegate et al. 2007;Greenwood 1992).This led Thorp and Covich (1991) to reject it from their index.The present results showed the seasonal stability of Tshred in our dataset, but analyses that are more specific could be performed to verify its stability regardless of the species considered, the development stage or even other ecological characteristics.
The only taxonomic composition metric selected, Abiv, is a long-standing bioindicator, dating back to the first international biomonitoring program (the Mussel Watch Program; (Be ´langer 2009).Bivalves are sensitive to disturbance and instability in marine ecosystems (Nerlovic ét al. 2011).They are additionally used as bioindicators to assess marine pollution along coastal zones (Yusof et al. 2004).The positive responses of 'Abiv' bivalves abundance to dry organic matter (DOM) as well as nutrients (Nitrate, TP, Ortho phosphates) are consistent, as bivalves are considered filterers of large volumes of water and also concentrators of organic matter (Petersen et al. 2019).However, the overabundance of organic matter can lead to a reduction in species abundance due to hypoxia and the accumulation of toxic by-products from the degradation of these materials (Hyland et al. 2005).This may explain the negative response of mussels to the high levels of organic pollution (BOD) and hypoxia (DO) observed.Ammoniacal nitrogen, as a by-product of TKN, is an important pollutant which is potentially harmful to aquatic organisms including bivalves and shellfish (Lu et al. 2022).Many bivalve species are known to inhabit muddy habitats (Lazo 2004) and also have a high sedimentation rate (Petersen et al. 2019).This explains the positive response of the ''Abiv'' metric for the ''Siltation'' pressure parameter.We can also explain the negative influence of dredging activities (Observation characteristic of this type of pressure proxy 'Depth') which likely lead to the suspension of sediment in the water column, by its abrasive effects on the eggs, by its reduction of the pumping rate of the bivalves and by their direct mortality (Wilber and Clarke 2001).In the context of Nokoue ´, resuspension is also driven by waves (wind effect or AWS) (Hounye `me `et al. 2023).The negative effects of acidification (Bressan et al. 2014;Ventura 2018), temperature and salinity on bivalves have also been demonstrated in the literature.Many studies have shown that bivalves react negatively to an increase in temperature above an optimal threshold (Bae et al. 2021;Galimany et al. 2011;Maynou et al. 2020;Steeves et al. 2018).It was the same for the salinity, which also beyond a certain optimal threshold would lead to mortality (Bae et al. 2021).The negative impact induced by the absence of bank macrophytes explained the disappearance of taxa living on their root systems.Low transparency (High-pressure level) corresponds to an increase in turbidity.The role of bivalves in increasing transparency is also noted (Rong et al. 2021).Mussel recruitment was negatively impacted by water color and turbidity (O ¨sterling and Ho ¨gberg 2014).As for the strong artificialization of soils (Oc_nN), its negative impact on bivalves could be explained in many ways.For example, agricultural pollution and pesticides can have negative effects on bivalves such as Corbicula fluminea (Dama ´sio et al. 2010).This becomes more obvious when we know that agriculture almost always has negative impacts and very few positive effects in marine environments (Galloway et al. 2008;OECD and Dı ´az 2010;Tilman et al. 2001).Last, we must note the systematic presence of bivalves such as corbulidae in the benthic macrofauna observed in the acadjas (Sankare et al. 2021).This would justify the positive response of 'Abiv' to the acadjas.
Selected diversity indices such as Simpson k index (Simpson 1949), Margalef index (Margalef 1958) and Berger-Parker index (Berger and Parker 1970), were already shown to be indicative of disturbance (Boyle et al. 1990;Ntislidou et al. 2018), supporting our results.These metrics, which have shown their usefulness in assessing the impact of water quality degradation (Teixeira 2010;Gamito et al. 2012;Nguyen et al. 2014;Van Loon et al. 2018;A ´vila et al. 2018) on phytoplankton (France ´et al. 2021) and fish, are also useful for invertebrates, according to our results.It should also be noted that water quality improvement is most often indicated by an increase in the Margalef index (Patrı ´cio et al. 2009), while the Berger-Parker index increases with the level of water pollution (Padmanabha 2011).The different responses of the metrics to the presence of acadjas seem to indicate a local increase in taxonomic and functional macrobenthic diversity or its heterogeneity.Note that all metrics respond consistently to one or another of the pressures, following the variabilities imposed on them according to the sensitivities of their constituent taxa.
Looking at the results, we note that among the relationships in the model in Fig. 4, the link between 'acadjas' and 'macrophytes' does not seem obvious.The fact that these two substrates are of the same nature, i.e. fixation supports for macroinvertebrates, possibly explains this causal relationship.
Literature or expert opinion is often used to define how parameters should respond to anthropogenic pressures (Drouineau et al. 2010(Drouineau et al. , 2012)).However, there are often contradictions between the meaning given by experts and the effects expected and induced by pressure proxies.The practice of extrapolation may explain this.The option taken in this study to consider only the valid inference of environmental effects produced by joint species distribution models (Roberts et al. 2022) seems appropriate given the results obtained here.The posterior probabilities of responses calculated from the available data provide trends, clear indications of the responses of each metric at each level of each pressure, and are valid in the absence of expert opinion.

Coherence of metrics reference values and stress thresholds
A comparison with recently proposed reference thresholds for West African rainforest rivers, based on best professional judgement (BPJ) and trisection methods (Kamelan et al. 2022), shows that our current stress thresholds are below their proposed lower limits for temperature (25.10-25.33°C) and nitrate (0.450-0.625 mg/L), and then place our stress thresholds above their upper limits for DO (6.33-6.40mg/L), pH (6.86-7.04)and TP (0.123-0.134 mg/L).Our approach leads to a more severe assessment, although we are not in the same environment.Bathymetric analyses of the bottom of the Nokoue ´lagoon, described since 1906 have positioned the depth around 2 m (Sossou-Agbo 2013).Considering that this depth value was already the result of alterations, a better compromise in relation to a stress threshold value would be a Depth [ 2 m.According to the triangle of division of the sedimentary dynamics, based on the method of Flemming (2000), it is obvious to consider the sandy texture (\ 5% mud) as the granulometric stress threshold value compared to the muddy one ([ 95% mud).The stress threshold value of 7.66% (Table 4) is of course close to the sandy texture.The threshold value of non-natural land use deducted from our analyse was high (45% for a minimum of 55% natural area).This is not very severe compared to the options taken in other studies/areas (Argillier et al. (2012).This may be due to the sampling points not being in a position to reflect the changes in the catchment (Hounye `me `et al. 2023).Regarding TKN, the 0.9 mg/L corresponds to the standard lower limit proposed for effluent discharge to surface waters, which is on the order of 100 mg/L (Dash 2013;Maiti 2003) and slightly below the minimum standards recommended by the US EPA.Knowing the anthropogenic nature of TKN, the proposed stress threshold value is in accordance with the values proposed by the US EPA standards.The stress threshold value obtained for total phosphorus (1.03 mg/L), it is well above the US EPA reference (0.03-0.1 mg/L).Such a high total phosphorus threshold is unique to this lagoon system which already has an endogenous stock of total phosphorus (Hounye `me `et al. 2023).Overall, the current stress thresholds observed are above the US EPA threshold for BOD (\ 6 mg/L for US EPA vs 7.19 mg/L) and orthophosphates (range of 0.025-0.18mg/L for US EPA vs 0.23-0.41mg/L).On the other hand, they are equivalent or better for nitrates (0.25-45 mg/L vs 0.25-0.41mg/L), DO (3-5 mg/L for US EPA vs 4.01-7.82mg/L), transparency (0.3-2.5 m for US EPA vs 0.56-1.32m), temperature (20-30 oC for US EPA vs 23.71-27.2) and pH (6-8.5 for US EPA vs 7.58-8.7).The proposed stress thresholds for DOM (2.75-4.18%)are in the same trend range as those of the ''good-moderate'' classes for the quality status obtained in the low salinity areas of the Mira estuary (Portugal) (Medeiros et al. 2012).For salinity in freshwater lagoons, the lowest stress threshold (2.08 %) was used as the stress threshold for the downstream freshwater part and the highest (8.86 %) for the upstream area close to the sea.These two stress thresholds appear to be consistent with the metrics that allowed their determination.Indeed, the Abiv metric, made up of freshwater taxa, corresponds to the lowest salinity stress threshold value, while the 'dMa' metric, which is very sensitive to natural variability, corresponds to the high stress threshold salinity value (Table 4).Supporting this argument, these two stress thresholds values fit into the top 2 of the five overlapping biological salinity zones (0-4, 2-14, 11-18, 16-27, and 24-marine).These salinity zones were defined through fish and invertebrate communities along mid-Atlantic estuaries (U.S.A.) (Bulger et al. 1993), then confirmed on estuaries in the northern Gulf of Mexico (Christensen et al. 1997), and finally on the Suwannee River (Florida) (Farrell et al. 2005).
Regarding the evidence on the metrics themselves, some consistent trends can be noted.Despite the increase in polychaetes (Serpulidae) under stressful conditions, they are still less sensitive to anthropogenic disturbance than crustaceans and molluscs (Reise 1982;Warwick and Clarke 1993;Wildsmith et al. 2009Wildsmith et al. , 2011)).Knowing this, the directions of quality and degradation shown by the 'Abiv' and 'Tperma' metrics at the 0.02 reference value seem appropriate.Some geographical areas with specificities can give reference values for the Berger-Parker index higher than 1 (Soudant and Belin 2011) as is the case (Table 4), whereas it should normally be in the range [0,1] (Patrı ´cio et al. 2009).The closer the Simpson index value k is to 0, the more diverse the community (Davis and Wing 2012).Thus, a more diverse ecosystem is indicated by the current reference value (\ 0.21) for the Sim_iL index.

Methodological considerations
The appeal to HMSC made in this approach is mainly justified by its outputs, the flexibility of its structure and its performance compared to other JSDMs (Ovaskainen et al. 2017).The approach developed has several advantages, including making the nature of interactions between stressors themselves and between stressors and their responses more explicit in the context of African ecosystems, where communities and their expected responses to stress are poorly documented.This would further improve our knowledge of these interactions (de Vries et al. 2021;Zhou et al. 2022;Brooks et al. 2021).Furthermore, the approach offers flexibility and potential for extension to many other biological systems, as it could incorporate expert advice, where available, in an explicit and reproducible manner, and also account for small sample sizes, high variability in transitional environments, and uncertainties in designing experimental protocols (Drouineau et al. 2010).The HMSC model precisely provides in its output an a posteriori probability for the different expected responses for the metrics, which once transformed are equivalent to expert knowledge.Depending on whether the credibility interval contains negative or positive values, it was also possible to assess the relevance and importance of the probability of a metric as a function of changes in the value of the environmental covariate causing the disturbance (Norberg et al. 2019).Thus, depending on whether the credibility interval contains negative or positive values, it can be said that a metric responds in a positive or negative way to the pressure parameter.Finally, the predictive power of the HMSC model through cross-validation is an important bias-reducing asset, as it allows for adjustment as the analysis progresses (Tikhonov et al. 2020).
The choice of the HMSC 3 (Step 2) model for the selection of the metrics is specific to this work, as the choice of the metrics best suited to each sampling period would have been sufficient.In fact, beyond the flexibility we were looking for in terms of seasonality, we felt that the manager should not burden his budget with extra sampling.As one of the assumptions in this work was to ensure that the sample points were independent, a random effect was introduced at the sample unit level (Ovaskainen et al. 2017).Our results clearly showed the negligible role of the sampling units in the model fit.Regarding potential generalization issues, the stratified random sampling option used in the HMSC 3 (step 2) model of metric selection seemed relevant to avoid pseudoreplication problems that could affect inferential validity (Lazic et al. 2020).Finally, the absence of autocorrelation (Fig. 3) is evidence that the independence assumption has not been violated.

Conclusions
The present results show that different stressors induce responses in the invertebrate communities sampled in the Nokoue ´freshwater lagoon.The approach used here is innovative because it allows us to dispense with expert opinion and to predict the responses of each of the metrics for any given stress level, which is essential for defining the reference values of the metrics and the stress thresholds.Our approach advocates the need to determine the sensitivity of the selected metrics to human perturbations as well as their robustness.The reference value and the derived stress thresholds are comparable to the existing ones, but with some specificities that could be attributed to the environment.The approach has proven its relevance in that it has produced very good inferences.The relevant metrics (Berger-Parker index, Margalef index, Simpson k index, % of substrate taxa preferring mud, % of taxa permanently attached to the substrate, % of shredded taxa, bivalve abundance), the reference values of the metrics and the stress thresholds provide at least a framework for assessing the ecological quality of this system with perspectives for improving management (pressure reduction measures, implementation of monitoring plans, restoration, prevention of erosion of its macrobenthic diversity, etc.).This assessment framework could be considered for extrapolation to other environments in the region, given the taxonomic composition of the selected metrics and the range extension of these species across the region.To the best of our knowledge, the set of response probabilities for each of the pressure levels indicating the trend of degradation, the reference values of the metrics and the stress thresholds obtained are particularly rare and therefore valuable data for African lakes or lagoons of this type.One of the prospects arising from the present work would be to envisage the construction of a multi-metric index and the valid extension of the approach to other ecosystems in the region.Indeed, with the selected metrics, we could envisage the construction of a multimetric index following the approach of Schoolmaster (2013), while initially defining disturbance gradients according to Hering et al. (2006).

Fig. 2
Fig. 2 Flowchart of the method for defining the reference thresholds (The dotted, rectangular or square representations are for the number of metrics selected and the circular ones specify the variant of the

Fig. 4
Fig. 4 DAG showing the relationships between environmental variables, stressors, forcing factors, and metrics: In yellow ellipse are the metrics, in rectangle the pressure proxies, in diamond shape the natural factors.The arrows interpret cause and effect

Fig. 5
Fig. 5 Positioning of the reference values (green) in relation to the distribution of the mean values of the metrics over the years (Red dotted line = reference value; blue dotted line = max value; orange dotted line = min value) Table 1 below shows the different pressure levels used to calculate the response probabilities of the different metrics.The intervals were defined on the basis of the change standards previously established for this lagoon (Hounye `me `et al. 2023).
By learnedBN (de Vries et al. 2021) and then extracting reference threshold values from probability density plots.

Table 1
Pressure levels and discretization ranges