Modeling drivers of biodiversity change emphasizes the need for multivariate assessments and rescaled targeting for management

The current policy and goals aimed to conserve biodiversity and manage biodiversity change are often formulated at the global scale. At smaller scales however, biodiversity change is more nuanced leading to a plethora of trends in different metrics of alpha diversity and temporal turnover. Therefore, large‐scale policy targets do not translate easily into local to regional management decisions for biodiversity. Using long‐term monitoring data from the Wadden Sea (Southern North Sea), joining structural equation models and general dissimilarity models enabled a better overview of the drivers of biodiversity change. Few commonalities emerged as birds, fish, macroinvertebrates, and phytoplankton differed in their response to certain drivers of change. These differences were additionally dependent upon the biodiversity aspect in question and which environmental data were recorded in each monitoring program. No single biodiversity metric or model sufficed to capture all ongoing change, which requires an explicitly multivariate approaches to biodiversity assessment in local ecosystem management.


| INTRODUCTION
Biodiversity change is a planetary concern. The 2019 global assessment by the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (IPBES), has raised awareness about the extent of the biodiversity crisis. The focus of the report is on species extinction, and this is also explicitly addressed by the Aichi targets of the Convention of Biological Diversity (Target 12, specifically;CBD, 2010). While the collective actions of society (and politics) failed to meet any of the 20 targets of Aichi at the global scale (CBD, 2020) the underlying assumption of a universal decline in biodiversity have been criticized as well Dornelas et al., 2014Dornelas et al., , 2018. The dynamics of biodiversity become much more nuanced at local to regional scales ; the scales at which most management decisions are made (Moloney et al., 2017). Different rates of immigration and extinction as well as shifting dominance result in a multifaceted local biodiversity change. This complexity of local biodiversity change cannot be channeled into simple univariate local biodiversity targets, which contrasts the global focus on rates of extinction. The fundamental challenges of local biodiversity management additionally comprise identifying relevant drivers of biodiversity change, as it is often these drivers (nutrient load, fragmentation, exploitation, and pollution) that actually can be targeted by management measures rather than biodiversity itself. While globally extinctions and threats can often be linked to certain causes (IPBES 2019), the long scientific history of illuminating biodiversity-environment relationships has taught us to expect different types of relationships. Single drivers such as eutrophication (Hillebrand et al., 2007) or warming Bastazini et al., 2021) triggered responses of different magnitude and even nature when comparing ecosystems and organism groups. These studies often focus on single or few drivers, but multivariate approaches show that diversity has different responses to different variables of environmental change (Grace et al., 2016). Whereas drivers of human impacts correlate rather strongly to each other in terrestrial systems, the same does not apply to marine ecosystems (Bowler et al., 2020). Thus, moving from global targets to local management options requires acknowledging the complexity of why biodiversity changes (multiple drivers and their interactions) and of how biodiversity changes (multiple facets of diversity and turnover). Consequently, monitoring data should be analyzed in a way that enables to document biodiversity change at the appropriate scales of target development and management in order to set conservation priorities accordingly.
Here, we use a well monitored ecosystem (Wadden Sea UNESCO World heritage site) as a showcase for additional inference that can be derived from combining different aspects of biodiversity and different potential drivers. Biodiversity change can either reflect the standing diversity (e.g., richness and evenness) of a locality (alpha diversity) or a region (gamma diversity) or the change in species composition between sites (spatial beta diversity) or time points (temporal turnover). Temporal trends in standing diversity can be disentangled from trends in temporal turnover (Dornelas et al., 2014, Hillebrand et al., 2018 and high compositional turnover can occur without strong trends in standing diversity if gains and losses balance each other over time (Brown et al., 2001). In contrast to the focus on the persistence or extinction of species at the global scale, local biodiversity change is additionally characterized by changes in the relative dominance of species (Hillebrand et al., 2008). While immigration and extinction rates affect the presence or absence of a species in a given place at a given time, environmental change also affects the relative abundance of species. Such altered dominance may even occur faster than actual extinctions or arrival of species and serve as a sentinel of biodiversity change.
For our study, we follow suggestions by Hillebrand et al. (2018) to analyze local to regional biodiversity change using a 2 Â 2 matrix, which combines measures of both standing diversity and temporal turnover that are or are not reflecting dominance change (Figure 1). For standing diversity, species richness (S) is the default measure if only absence and presence of species shall be reflected. Although simulations clearly show that S is easily affected by systematic errors induced by comparing different sampling efforts, species pool sizes, and spatial distribution of individuals (Chase & Knight, 2013). From the plethora of dominance-weighted indices of standing diversity, the effective number of species (ENS) turned out to be relatively robust against sampling and statistical issues (Chase & Knight, 2013). ENS is the number of equally abundant species that would be needed to achieve the same diversity as observed. If in a local assemblage, all species are equally abundant, then ENS equals S, in all other cases ENS < S.
For turnover, metrics capture the relative difference in composition. Jaccard dissimilarity is a presenceabsence based metric ranging from 0 to 1 (0 = no difference in species composition and 1 = no overlap). As S for alpha diversity, Jaccard gives equal weight to rare and abundant species. For dominance weighted turnover, again a plethora of indicies has been derived. We opted for the less well-known Wishart dissimilarity (Wishart, 1969) as it emphasizes dominant species to the same degree as ENS does. Thus, we represent temporal turnover by two species exchange ratios, one based on richness (SERr; i.e., presence-absence, SERr = Jaccard Dissimilarity) and one based on abundance (SERa = Wishart's dissimilarity). However, it should be noted that this article is not about advocating a certain metric but to elucidate how much more we can learn about the facets of local biodiversity change when using the 2 Â 2 matrix of presenceabsence and dominance driven indices of alpha diversity and turnover.
In this context, we also focus on the identification of potential drivers for biodiversity change in a multivariate context, as the temporal trends per se have been reported elsewhere. Rishworth et al. (2020) demonstrated how using these different biodiversity measures captured aspects of temporal biodiversity change for multiple organism groups in two study sites, the German Wadden Sea, and the South African Algoa Bay. Here, we analyze the role that several environmental variables play in biodiversity change in the Wadden Sea by combining two statistical modeling approaches: Structural equation modeling (SEM) for analyzing alpha diversity, generalized dissimilarity modeling (GDM) for temporal turnover. Our analysis relies on data from standard monitoring of fish, birds, macrozoobenthos, and phytoplankton from Germany and The Netherlands. Thereby, we aim at answering the following questions: 1. Can biodiversity change be described by single aspects of standing diversity or temporal turnover, or does it only fully emerge from a multifaceted approach? 2. Do different environmental drivers affect different biodiversity metrics similarly and consistently across organism groups or are these relationships between driver and response specific?
The Wadden Sea, as a region of the North Sea, stretches along the coastlines of the Netherlands, Germany, and Denmark and is the world's largest intertidal mudflat system ( Figure S1). We analyzed community composition data by organism group from various monitoring programs of the Wadden Sea ( Figure 1). We drew the data from established local scale surveys of trawled fish, migratory birds, benthic macroinvertebrates (hereafter macrozoobenthos), and phytoplankton, each reporting the presence and abundance of species (rarely higher taxonomic levels) in the Wadden Sea of Germany (specifically Lower Saxony) and the Netherlands (Table 1). Fish were caught at several stations ( Figure S1) with beam trawls, mainly in March and July (Wetjen et al., 2014). Migratory birds were counted fortnightly or monthly in permanent counting area covering almost the entire coastline ( Figure S1; Kleefstra et al., 2022). Phytoplankton sampling was done annually at least during the vegetation period (March-September/October) at representative monitoring sites ( Figure S1) and at least once a month (NLWKN, 2013). Macrozoobenthos was sampled at fixed stations by extracting a sediment sample at least every 6 months (Drent et al., 2017). To ensure consistency per organism group across countries, we checked the data for consistent naming and cleaned accordingly. The datasets were in themselves consistent in sampling, laboratory analysis, and taxonomic identification. To focus on long-term annual changes in biodiversity patterns over time, we calculated diversity measures first F I G U R E 1 Flow diagram of our general approach to biodiversity models of four Wadden Sea organism groups: Fish, birds, macrozoobenthos, and phytoplankton and then took their yearly averages of every site for each of the datasets. We then matched the resulting biodiversity data with yearly environmental data averages collected at the same site and over the same time frame. This generally resulted in an averaging out of extreme events, such as heatwaves or extremely cold winters. Sometimes, the environmental or species data presented considerable gaps, which were worsened when matching both sets together. The resulting sample size after eliminating the gaps varied for each of the four organism datasets (Table 1). More analysis details are in the Supporting Information).

| Statistics
The two levels of inference, standing diversity and turnover, required two different statistical approaches ( Figure 1). For standing diversity, the annual mean S and ENS were related to the annual mean in environmental variables as well as annual average biomass. For temporal turnover, the analysis needs to be performed in a dissimilarity perspective, that is, the change in composition related to the gradient in the environmental variables. It should be noted that we focused on temporal turnover between adjacent years, that is, altered as SERr and SERa with the same delta-time of 1 year. For analyzes of cumulative turnover and species composition see Rishworth et al. (2020). In both parts, we used two measures for biodiversity and for species turnover, one based on presence/absence data and the other on relative species abundance in samples ( Figure 1). For standing diversity, we used S calculated as the number of species in a sample and the ENS calculated as the inverse Simpson index (Hillebrand et al., 2018). We conducted all analyzes in R version 4.0.3 (R-Core Team, 2021). We calculated both S and ENS using the vegan-package (Oksanen et al., 2020). The turnover indices SERr and SERa were calculated with function turnover: https://github.com/ AlexRyabov/turnover. All scripts and corresponding data are available on our GitHub-repository: https://github. com/JanDajka/BiodivModelsWaddenSea.

| Standing biodiversity models
We analyzed the effects of environmental variables on standing diversity separately for each of the organism groups using piecewise SEMs implemented in the piecewise SEM-package (Lefcheck, 2016). A SEM weaves multiple linear relationships together, allowing us to find a set of interacting response variables and to display the network of direct and indirect causal links between them (Grace, 2006). The advantage of piecewise SEMs over traditional SEMs is that we can add random effects, nested structures, and non-normal distributions to the models (Lefcheck & Duffy, 2015). Classically known "response variables" are here "endogenous variables" and environmental predictors are often called "exogenous variables." We fitted linear models with generalized least squares (GLS) using the nlme-package (Pinheiro et al., 2017) to our data and used latitude and longitude of the sites nested in year as correlation structures to account for spatial-temporal autocorrelation. We used the following overall GLS-structure: For each of the four organism groups' SEMs, we created GLS-models for S, ENS, and log-transformed biomass. We scaled all explanatory variables to a mean of 0 and standard deviation of 1 to allow for a comparisons of T A B L E 1 Data used in final model construction for all four organism groups of the Wadden Sea; institutions: Lower Saxon Wadden Sea National Park Authority (NLPVW), Royal Netherlands Institute for sea research (NIOZ) lower Saxon State Department for waterway, coastal and nature conservation (NLWKN), directorate-general for public works, and water management (Rijkswaterstaat) effect sizes (Schielzeth, 2010). Our model validation showed a normal distribution of the residuals and no clear patterns between residuals and fitted values for all linear models, indicating that model fit could not be improved further. The resulting linear models were then woven together into our final piecewise SEM. Variables that were strongly correlated (abs[r] > 0.5) were declared as correlated in the SEM-structure. In order to estimate the goodness of fit of the piecewise SEM with a Fisher's C-test, we removed the least significant variable of each linear model as is commonly done for SEMs (Maureaud et al., 2019). The above-described process was the same for all organism groups, just with different environmental variables (Table 1) and without biomass for birds since this dataset is based on counting resting birds rather than collecting and weighing them.

| Temporal turnover models
We analyzed the dependence of the annual species turnover on the gradient in environmental variables using GDM with the gdm-package (Manion et al., 2017). We decided against also applying SEM to turnover due to the nature of turnover data being dissimilarities between years. The resulting graphs would have shown the contribution of differences in environmental variables to turnover and made interpreting results difficult. Instead, GDM relates changes in temporal turnover to gradients in environmental variables using monotonic I-spline functions (Rillo et al., 2022). In the resulting plots, the maximum height of each curve displays how much overall turnover is explained by each environmental variable. A leveling-off of the curve shows where along the environmental gradient most (or least) of change in turnover occurs. GDM is typically used to explain the effect of environmental variables on spatial differences of species diversity (β-diversity; Ferrier et al. 2007;Latombe et al., 2017). We applied GDM to temporal differences of species diversity by using a matrix showing year to year turnover in species identities (SERr) and species composition (SERa) for every site for which we had mean annual measurements of environmental variables. All environmental variables were scaled variables to a mean of 0 and standard deviation of 1 to allow direct comparison between sites.

| RESULTS
For each organism group, SEM and GDM analyzes allow us to identify a characteristic set of environmental drivers of diversity. However, the strength and nature of each environmental factor's influence strongly depend on the organism group and the biodiversity metric in question. Each of the highlighted relationships can be analyzed in depth individually, which is useful for management if the impact of a certain environmental variable on biodiversity needs to be checked. Because this would make our results section incredibly dense, we focus only on the main differences highlighted by modeling the 2 Â 2 matrix. In fish, the key effects of environmental variables highlighted were quite comparable between metrics (Figure 2). Increasing temperature had a positive effect on both measures of standing diversity, ENS and S, but no direct effect on biomass. S and ENS were negatively correlated. Instead, biomass was positively associated to increasing S and negative to increasing ENS. The temperature gradient also was the strongest predictor of dissimilarity, whereas pH and oxygen had weak impacts in both SEM and GDM.
For birds, the SEM did not detect any effects of environmental variables on S or ENS, which however were positively correlated (Figure 3). In contrast, the GDM revealed that temporal dissimilarity in bird composition increased with larger gradients in temperature, precipitation, and macrozoobenthos biomass that served as a proxy of food availability.
Macrozoobenthos showed the same positive relationship between S and biomass as fish and a weak positive correlation between S and ENS. The SEM and GDM showed differential impact patterns of environmental variables on the biodiversity metrics ( Figure 4). The SEM showed higher explained variance in S than in ENS. Positive responses of S to grain size and salinity but negative to temperature were only partially matched in ENS, which increased with both salinity and temperature, but was unaffected by grain size. By contrast, the GDM revealed stronger response in the dominance weighted SERa than in SERr. However, in both cases grain size showed the strongest impact, salinity the least.
For phytoplankton, the SEM explained slightly more variance in S than ENS whereas the GDM showed clearer patterns for SERa than SERr ( Figure 5). Moreover, standing diversity strongly responded to silicon and temperature, whereas nutrients and pH played a stronger role in the GDM. S decreased with silicon while ENS increased with silicon while both S and ENS increased with temperature. Increasing nitrogen had a very strong effect on ENS and none on S. Phytoplankton biomass increased with phosphorus and pH. Increasing silicon, temperature, and salinity negatively affected biomass. Moreover, increasing S positively affected biomass while ENS had a strong negative effect on biomass.

| DISCUSSION
Biodiversity change is nuanced, especially at local scales . Our study confirms that this also applies to the environmental variables driving the changes in biodiversity. First, our study emphasizes the importance of using multivariate assessments of biodiversity. Simple univariate assessments would miss a large part of both biodiversity change itself and the impact of different environmental variables. This also has important ramifications for the ways in which the targeting efforts of management function. Effective management relies on effective targets-from which to measure change and meet goals. These in turn need to be established from data that enable understanding of biodiversity change and be reflective of their drivers as well. As our findings show, multivariate assessments at local scales give a lot more insight into the complex nature of F I G U R E 2 Effects of environmental factors on diversity metrics for fish of the Wadden Sea; estimates for standing diversity and biomass are modeled by SEM (a): Endogenous variables with arrows going into them; exogenous variables with arrows coming out of them; dashed arrows: insignificant relationship (p > .1); solid arrows: significant relationship; thicker arrows: higher significance; red: negative relationship; black: positive relationship; gray: indiscriminate relationship (coefficient value below 0.1 or À 0.1); numbers: coefficient strength (from À1 to 1, strengths of individual arrows can be compared with each other); dotted connectors without arrows: correlations (with correlation coefficient numbers) between variables; environmental effects on temporal turnover estimated with GDMs (b) richnessbased SERr (top) and abundance-based SERa (bottom): Each line is a significant site, sites where certain links to variables were insignificant, had a coefficient of <0.1 or the deviance explained by the model was <15%, were not plotted; environmental variables are scaled to a mean of 0 and standard deviation of 1 to make their impacts directly comparable, the height of the lines signifies a higher impact on temporal turnover of the respective organism community, leveling-off of the curve shows where along the environmental gradient most (or least) of change in turnover occurs biodiversity change (Rishworth et al., 2020). Additionally, our study shows how multivariate measures of biodiversity change (Hillebrand et al., 2018) can be linked to environmental variables using SEM and GDM. The primary strength of both modeling approaches is that the impacts of all modeled variables are directly comparable and conservation priorities as well as suitable conservation measures can be formulated along the links highlighted in our study. Through these links, we move closer to enabling management of biodiversity change at local scales, where tailored approaches will be more effective. At these smaller scales, targets may be more realistic and meaningful to set and reach.
Our assessment approach requires recording species identities, their abundance, their biomass, and respective environmental variables. We want to urge management to improve their efforts in environmental monitoring to not underestimate the extent of biodiversity change.
F I G U R E 3 Effects of environmental factors on diversity metrics for birds of the Wadden Sea; estimates for standing diversity and biomass are modeled by a SEM (a): Endogenous variables with arrows going into them; exogenous variables with arrows coming out of them; dashed arrows: insignificant relationship (p > .1); solid arrows: significant relationship; thicker arrows: higher significance; red: negative relationship; black: positive relationship; gray: indiscriminate relationship (coefficient value below 0.1 or À0.1); numbers: coefficient strength (from À1 to 1, strengths of individual arrows can be compared with each other); dotted connectors without arrows, correlations (with correlation coefficient numbers) between variables; environmental effects on temporal turnover estimated with GDM (b) richness-based SERr (top) and abundance-based SERa (bottom): each line is a significant site, sites where certain links to variables were insignificant, had a coefficient of <0.1 or the deviance explained by the model was <15%, were not plotted; environmental variables are scaled to a mean of 0 and standard deviation of 1 to make their impacts directly comparable, the height of the lines signifies a higher impact on temporal turnover of the respective organism community, leveling-off of the curve shows where along the environmental gradient most (or least) of change in turnover occurs A clear limitation in our study, for instance in the SEM for birds, arose from the lack of available data, especially those of environmental variables. Currently, the level of detail and sampling consistency of most Wadden Sea species data (especially from Lower Saxony) far outmatch those of corresponding environmental variables. Our assessment here provides a method on how to model the question of biodiversity change and its drivers but does not provide a final answer to it. To continually observe change, monitoring concepts need to consistently match more extensive environmental data to species data and, ideally, collect them together. F I G U R E 4 Effects of environmental factors on diversity metrics for macrozoobenthos of the Wadden Sea; estimates for standing diversity and biomass are modeled by SEM (a): Endogenous variables with arrows going into them; exogenous variables with arrows coming out of them; dashed arrows: insignificant relationship (p > .1); solid arrows: significant relationship; thicker arrows: higher significance; red: negative relationship; black: positive relationship; gray: indiscriminate relationship (coefficient value below 0.1 or À0.1); numbers: coefficient strength (from -1 to 1, strengths of individual arrows can be compared with each other); dotted connectors without arrows: correlations (with correlation coefficient numbers) between variables; environmental effects on temporal turnover estimated with GDM (b) richness-based SERr (top) and abundance-based SERa (bottom): Each line is a significant site, sites where certain links to variables were insignificant, had a coefficient of <0.1 or the deviance explained by the model was <15%, were not plotted; environmental variables are scaled to a mean of 0 and standard deviation of 1 to make their impacts directly comparable, the height of the lines signifies a higher impact on temporal turnover of the respective organism community, leveling-off of the curve shows where along the environmental gradient most (or least) of change in turnover occurs F I G U R E 5 Legend on next page.
Better protection and management concepts can indeed be achieved through improved biodiversity monitoring, data analyzes and modeling approaches. Approaches like ours can be part of the blueprint for these concepts. With increasing sample size from longer time series and from fitting species data with more extensive environmental data, modeling approaches such as SEM and GDM become more accurate, can detect more links, and thus provide a better picture of the modeled system. The phytoplankton dataset in our analysis was paired with more extensive environmental data and we were thus able to create a more detailed model compared to the other organism groups. In contrast, macrozoobenthos species data were measured with higher temporal and spatial resolution but was only complemented with three environmental variables that were congruent between Germany and the Netherlands. This was the also the case for the bird data, except that a congruent species-environmental-data combination was only available for Germany. The fish data presented themselves with the lowest sample size from Germany only, which is likely the reason for the few significant relationships in the models. The Wadden Sea stretches across three countries and more extensive datasets are especially hard to come across due to its transnational geography and trilateral governance. Cross-national agreements on environmentally matched, continuously collected time series data will enhance the complexity and accuracy of models in the future, and in turn effective management regimes.
The implications of biodiversity change are being increasingly well understood, but, due to their complexity, cannot be easily addressed. Recently proposed efforts to pursue univariate, global, targets for biodiversity comparable to the 2.0 C-climate target (Rounsevell et al., 2020) would lead to ineffective governance of this complex problem. It is evident from our study and the wider literature (Hillebrand et al., 2018;Rishworth et al., 2020) that multivariate biodiversity assessments at local scales-the scales that are vital for management decisions-should be pursued instead. Focusing on only singular biodiversity metrics in one-dimensional targets will miss major developments of biodiversity change, multi-dimensional targets must be used instead. These multidimensional targets can be based on multivariate assessments such as ours. The four biodiversity metrics S, ENS, SERa, and SERr are based on species identity and abundance-data that are often already being collected in many monitoring programs. We added biomass to our assessments, but this is not a necessity. For practitioners, this means that multivariate assessments are then only an additional calculating effort. Instead of picking out a single biodiversity metric to focus on, biodiversity monitoring needs to consider at least these four and ideally monitor environmental variables alongside it. Next steps could include connecting model results such as ours to nature's contributions to people (Kadykalo et al., 2019) to reveal the direct effect of biodiversity change on human communities. This provides a chain of links that policy makers can follow to effectively align their local scale management priorities along meaningful multivariate assessments of biodiversity.

AUTHOR CONTRIBUTIONS
Jan-Claas Dajka conceived core ideas and design of the study, harmonized the data, analyzed the data, and wrote the first draft of the manuscript. Josie Antonucci di Carvalho harmonized the data interpreted the results. Gregor Scheiffarth made the data available. Lena Rönn made the data available. Rob Dekker made the data available. Bo Leberecht interpreted the results. Helmut Hillebrand conceived core ideas and design of the study and interpreted the results. Alexey Ryabov analyzed the data and interpreted the results. All authors contributed critically to subsequent drafts and gave final approval for publication.

ACKNOWLEDGMENTS
We are grateful for the data provision by Nationalparkverwaltung Niedersächsisches Wattenmeer (NLPVW) and Niedersächsischer Landesbetrieb für Wasserwirtschaft, F I G U R E 5 Effects of environmental factors on diversity metrics for phytoplankton of the Wadden Sea; estimates for standing diversity and biomass are modeled by SEM (a): Endogenous variables with arrows going into them; exogenous variables with arrows coming out of them; dashed arrows: insignificant relationship (p > .1), solid arrows: significant relationship; thicker arrows: higher significance; red: negative relationship; black: positive relationship; gray: indiscriminate relationship (coefficient value below 0.1 or À 0.1); numbers: coefficient strength (from À1 to 1, strengths of individual arrows can be compared with each other); dotted connectors without arrows: correlations (with correlation coefficient numbers) between variables; environmental effects on temporal turnover estimated with GDM (b) richness-based SERr (top) and abundance-based SERa (bottom): Each line is a significant site, sites where certain links to variables were insignificant, had a coefficient of <0.1 or the deviance explained by the model was <15%, were not plotted; environmental variables are scaled to a mean of 0 and standard deviation of 1 to make their impacts directly comparable, the height of the lines signifies a higher impact on temporal turnover of the respective organism community, leveling-off of the curve shows where along the environmental gradient most (or least) of change in turnover occurs