Influence of land-derived stressors and environmental variability on compositional turnover and diversity of estuarine benthic communities

It can be challenging to differentiate community changes caused by human activities from the influence of natural background variability. Using gradient forest analysis, we explored the relative importance of environmental factors, operating across multiple spatio-temporal scales, in influencing patterns of compositional turnover in estuarine benthic macroinvertebrate communities across New Zealand. Both land-derived stressors (represented by sediment mud content and total sediment nitrogen and phosphorus content) and natural environmental variables (represented by sea surface temperature, Southern Oscillation Index, and wind−wave exposure) were important predictors of compositional turnover, reflecting a matrix of processes interacting across space and time. Generalized linear models were used to determine whether measures of benthic macroinvertebrate diversity, which are commonly used as indicators of ecological health on a local scale, changed in a way that was consistent with the compositional turnover along the environmental gradient. As expected, compositional turnover along landderived stressor gradients was negatively associated with diversity indices, suggesting a decline in ecological health as land-derived stressors increase. This study moves towards an ecosystem-based management approach that focusses on cumulative effects rather than single stressors by considering how multiple land-derived stressors influence indicators of estuarine health, against a background of natural variability across several spatio-temporal scales. OPEN ACCESS Benthic macroinvertebrate communities can indicate the health of an estuary and are influenced by both landderived stressors and natural variability. Photos: Cawthron Institute


INTRODUCTION
Understanding the influence of human activities on coastal ecosystems requires the separation of natural and anthropogenic sources of environmental variability. Partitioning these effects is particularly diffi-ABSTRACT: It can be challenging to differentiate community changes caused by human activities from the influence of natural background variability. Using gradient forest analysis, we explored the relative importance of environmental factors, operating across multiple spatio-temporal scales, in influencing patterns of compositional turnover in estuarine benthic macroinvertebrate communities across New Zealand. Both land-derived stressors (represented by sediment mud content and total sediment nitrogen and phosphorus content) and natural environmental variables (represented by sea surface temperature, Southern Oscillation Index, and wind−wave exposure) were important predictors of compositional turnover, reflecting a matrix of processes interacting across space and time. Generalized linear models were used to determine whether measures of benthic macroinvertebrate diversity, which are commonly used as indicators of ecological health on a local scale, changed in a way that was consistent with the compositional turnover along the environmental gradient. As expected, compositional turnover along landderived stressor gradients was negatively associated with diversity indices, suggesting a decline in ecological health as land-derived stressors increase. This study moves towards an ecosystem-based management approach that focusses on cumulative effects rather than single stressors by considering how multiple land-derived stressors influence indicators of estuarine health, against a background of natural variability across several spatio-temporal scales. cult in estuaries, due to the inherent complexity of these ecosystems, which are highly variable in both space and time (Elliott & Quintino 2007, Dauvin & Ruellet 2009). The impact of human activities is often assessed using benthic macroinvertebrate communities because they cover numerous trophic levels, exhibit different stress-tolerances, and can integrate the effects of multiple stressors over time (Pearson & Rosenberg 1978, Dauer 1993, Borja et al. 2000. These animals are also an important component of estuarine systems, playing essential roles in ecosystem structure and function (e.g. nutrient cycling, energy transfer to higher trophic levels, sediment stabilization; Snelgrove 1997, Levin et al. 2001, Lohrer et al. 2004). However, it can be challenging, particularly at large scales, to differentiate community changes caused by stressors from the influence of strong natural environmental gradients.
Estuarine benthic community structure is influenced by a range of natural, temporally varying factors, that operate at local (e.g. wind−wave exposure, sediment grain size, salinity, and predation; Snelgrove 2001) and broad (e.g. temperature, climate patterns; Engle & Summers 1999, Hewitt et al. 2016, Denis-Roy et al. 2020) spatial scales. Many of these natural factors can also be considered anthropogenic stressors when they exceed their natural range of variation as a result of human activities (Sanderson et al. 2002, Halpern et al. 2007). Estuarine communities are often exposed to multiple and cumulative stressors, and these commonly interact in multiplicative and non-linear ways (Crain et al. 2008, Darling & Côté 2008, deYoung et al. 2008). Many of these stressors are diffuse, operating in incremental stages and often over broad scales, particularly land-derived stressors like sedimentation and nutrient loading.
Land-derived stressors often represent natural processes that have been greatly accelerated by human activities and begin to have negative effects when their rate of delivery exceeds the assimilative capacity of the system. Sedimentation and nutrient loading are recognized as major threats to the health and functioning of estuaries globally (NRC 2000, EU Marine Strategy Framework Directive 2008, Magris & Ban 2019. For example, sedimentation rates have increased by 1 to 2 orders of magnitude in some places (Thrush et al. 2004), and 30−60% of estuaries in the USA and Europe are affected by nutrient enrichment (Bricker et al. 2008, EEA 2012. Adverse effects arising from these stressors (e.g. smothering of benthic communities, reduction in water quality; Valiela et al. 1992, Ellis et al. 2002 often manifest as reductions in species richness, evenness, and diver-sity, and a loss of rare taxa (Smith & Kukert 1996, Lardicci et al. 1997, Tagliapietra et al. 1998, Thrush et al. 2003b, Ellis et al. 2004. Community changes caused by land-derived stressors in the short term are often subtle, but cumulatively, across large spatio-temporal scales, they can drive substantial disruptions to ecosystem functioning. Sometimes stressors can also cause abrupt shifts in ecosystem functioning if a tipping point is reached (Hewitt & Thrush 2019) or in response to extreme pulse disturbances (e.g. storms; Thrush et al. 2003a).
Disentangling this complex web of factors is critical for understanding whether observed changes in benthic communities are indicative of degradation in ecosystem health or merely a result of natural environmental variation. Although the need to account for natural variability has been identified in regulatory documents such as the EU Water Framework Directive (2000), and is integral to ecosystem-based management (Arkema et al. 2006), it is seldom incorporated into assessment protocols due to the difficulty of teasing these factors apart and a perceived need to keep things simple (Irvine 2004). The influence of stressors and environmental variables operating on local scales needs to be considered within the context of processes acting across broader geographic and time scales within which the community is embedded (Ricklefs 1987). Such studies are uncommon in estuaries (but see Hewitt & Thrush 2009, de Juan & Hewitt 2011, Denis-Roy et al. 2020 because they require good spatio-temporal data along with methods that can quantify community response across multiple environmental gradients, while accounting for potential non-linearity and interactions among environmental variables. New Zealand spans 3 water masses, 15 degrees of latitude, and a variety of estuary types, providing an ideal place to investigate estuarine community responses under a range of environmental conditions. We used gradient forest (GF) analysis  to separate natural and anthropogenic drivers of benthic macroinvertebrate compositional turnover using a large nation-wide estuary monitoring dataset. In particular, we were interested in the community response to 2 pervasive land-derived stressors acting at a local (site) scale (sedimentation and nutrient loading) and 3 natural environmental variables representing both broad-scale (national) climate fluctuations (sea surface temperature [SST] and Southern Oscillation Index [SOI]) and local-scale processes (wind−wave exposure). Although we have classified these variables as either land-derived stressors or natural environmental variables for the purposes of this study, we acknowledge that they could be considered as either natural components of the system or human-induced stressors, depending on values relative to background levels.
GF is one of several of statistical approaches that can be used to model constrained relationships between communities and their environments (e.g. canonical correspondence analysis, multivariate regression trees, generalized dissimilarity modelling; reviewed by Ferrier & Guisan 2006). It has been used to explore marine ecosystem responses to anthropogenic and environmental pressures (e.g. Large et al. 2015, Samhouri et al. 2017, Couce et al. 2020 because it can model non-linear response shapes, deal with correlated predictors, and incorporate complex interactions between multiple predictors . It does this by combining information from multiple tree-based regression models (random forests, RF), one for each taxon, to provide a measure of compositional turnover across environmental gradients. Compositional turnover, sometimes referred to as beta diversity, is the component of regional biodiversity that accumulates due to inter-site variation in local species assemblages (Anderson et al. 2011, Socolar et al. 2016. Examining patterns in compositional turnover is important for identifying and un derstanding the processes that maintain species diversity across large spatial and temporal scales (Ricklefs 1987, Soininen 2010. For example, the large natural environmental gradients in SST and wind−wave exposure across our New Zealand-wide dataset would be expected to generate changes in turnover as community composition changes on a local scale. Although compositional turnover provides us with a measure of change in benthic communities in response to different environmental variables, it does not provide information on whether these changes translate into positive or negative effects on benthic communities on a local scale (Socolar et al. 2016). For example, the early stages of anthropogenic impact may cause localized species losses, leading to an increase in compositional turnover. However, anthropogenic impacts can also reduce compositional turnover rates, such as occurs when bottom-trawling destroys microhabitats leading to homogenization of benthic communities (Hewitt et al. 2005). Therefore, we used generalized linear models (GLMs) to determine whether measures of benthic macroinvertebrate diversity (i.e. species richness, evenness, diversity, and numbers of rare taxa), which are commonly used as indicators of ecological health on a local scale, changed in a way that was consistent with the compositional turnover along the environmental gradient.
We hypothesized that: (1) Both land-derived stressors and natural environmental variables will be important in predicting patterns of compositional turnover in estuarine benthic macroinvertebrate communities, reflecting a matrix of processes acting at different scales; (2) Compositional turnover along land-derived stressor gradients will have a negative relationship with species richness, evenness, diversity, and numbers of rare taxa.

Study sites
Data were obtained from estuarine monitoring surveys undertaken between 2001 and 2017 by New Zealand's regional government authorities (334 site/times sampled across 208 sites in 34 estuaries; Berthelsen et al. 2020a,b). The study sites ( Fig. 1; and see Table S1 in the Supplement at www. int-res. com/ articles/ suppl/ m666 p001 _ supp. pdf) spanned 12 degrees of latitude, 3 geomorphological estuary types (tidal lagoons, shallow drowned valleys, deep drowned valleys; Hume et al. 2016), and a wide spectrum of land-use intensities. Samples were collected between November and May (late spring to autumn), with the majority (70%) collected during the austral summer. Surveys were generally carried out according to a standardized protocol (Robertson et al. 2002), with samples collected from sites located at mid-tolow tidal height away from point-source discharges. To standardize for salinity effects, sites suspected to be significantly influenced by freshwater, based on proximity and flow rate of nearby streams, were removed from the dataset as well as any sites located within freshwater-dominated estuaries (i.e. tidal river mouth estuaries; Hume et al. 2016).

Macroinvertebrate data
Benthic macroinvertebrate samples (n = 3−15 replicates per site/time) were collected using a 13 cm diameter core, extending 15 cm into the sediment, and sieved to 500 μm. Experts identified organisms to the lowest practicable resolution. Taxonomic nomenclature followed the World Register of Marine Species (WoRMS Editorial Board 2017), and where differences in taxonomic resolution arose, we aggregated to higher taxonomic groups. This taxonomic aggregation may have obscured some of the true diversity; however, as taxa from all sites/times were treated the same, diversity indices are comparable on a relative scale. Some taxa were removed from the dataset before analysis, including taxa not wellrepresented by this sampling method (e.g. Bryozoa, meiofaunal taxa), those identified to relatively coarse taxonomic groups (e.g. Polychaeta, Annelida), larval planktonic groups (e.g. megalopes, eggs), non-marine taxa (e.g. Insecta, Acari), vertebrates, plants, and bacteria. Most (74%) of the remaining 122 taxa were identified to genus or species level. Abundance data were used in all analyses, with data from replicate macroinvertebrate samples averaged by site/time.

Environmental variables
Data for environmental variables (land-derived stressors and natural environmental variables) considered potentially important for influencing estuarine benthic macroinvertebrate turnover were collected concurrently with macroinvertebrate samples or collated from existing datasets ( Table 1). As community responses reflect environmental processes operating over a range of scales (Thrush et al. 2005), these variables were chosen to incorporate both local-scale factors that varied by site and broad-scale factors that operated at an estuary or national scale. We limited our assessment to 6 environmental variables, as the inclusion of many variables in regression tree approaches (such as GF) has been shown to provide only minimal improvement in predictive accuracy and to complicate interpretation of model outcomes (Leathwick et al. 2006).

South Island
North Island percentage of < 63 μm out of the < 2 mm sediment fraction) because the maximum grain size differed between analysis methods (e.g. Malvern Mastersizer laser only analyses grains < 2 mm, while all grain sizes are usually analysed during wet sieving). Mud concentrations at site/times used in this study covered the full spectrum (0−99% mud content), with a median of 13% (Table 1). Sediment total nitrogen (TN) and total phosphorus (TP) were used as proxies for nutrient loading. Despite slight variations in methods used to analyse TN and TP at different sites, results were assumed to be generally comparable by Berthelsen et al. (2020a). Values less than the analytical detection limit (ADL) were assigned values of ADL/2. Nutrient concentrations provided a wide stressor gradient, with maximum TN (4133 mg kg −1 ) and TP (1836 mg kg −1 ) values comparable to highly polluted estuarine sites worldwide (Oviatt et al. 1984, Gillespie & MacKenzie 1990, Sánchez-Moyano et al. 2010, Cao et al. 2011; Table 1). However, higher sediment TN values have been observed in some European estuaries (e.g. up to 8600 mg kg −1 in Bilbao Estuary, Spain; Saiz-Salinas 1997), and median values for both these nutrients were relatively low across the study sites (410 mg kg −1 TN, 340 mg kg −1 TP). Data from replicate sediment samples (mud, TN, and TP) were averaged by site/time.

Stewart Island
SST, SOI, and wind−wave exposure are natural environmental variables known to influence estuarine biodiversity (Engle & Summers 1999, Hewitt & Thrush 2009, Hewitt et al. 2016, Denis-Roy et al. 2020. SOI is a measure of the strength of the El Niño−Southern Oscillation, which occurs every 2 to 7 yr and is the largest source of natural variability in the global climate (Diaz 2005). Monthly estimates of SOI, corresponding with each site/time, were used as a measure of broad-scale temporal variability in climate. While extreme values in our 16 yr dataset were slightly less than those observed over more extended periods (−3.6 to 3.3 range since 1882), the dataset captured both El Niño and La Niña events (prolonged monthly average SOI below −1 or above 1, respectively).
Modelled average monthly SST data were obtained from the Jet Propulsion Laboratory Multiscale Ultra-high Resolution Sea Surface Temperature Project (NASA/JPL 2015) as a broad-scale measure of temporal and spatial variability across the study area. Values were taken from a location near the seaward entrance of each estuary and corresponded with the month and year of macroinvertebrate and sediment sample collection. Where SST data were not available for a site/time (n = 23), median SST across other site/times and within the same estuary, or a nearby estuary, was used. Median values were calculated from site/times sampled in the same month as the missing site/time SST where possible.
Wind−wave exposure was calculated for each site following a topographical method similar to that developed by Burrows et al. (2008). Wind direction and speed data, across 3 years of records, were obtained from the nearest regional airport and predominant winds binned into 45° intervals to give a measure of wind−wave disturbance from 8 directions. Around each site, the distance to land (fetch, measured in m) was calculated for every 1°, and each fetch value was multiplied by the total number of days when the predominant wind was from that direction and the wind speed (surface wind at 09:00 h, m s −1 ) for those days. Outputs were divided by 100 000 to convert the data to a smaller scale. Where sites were too close to land to calculate expo- sure metrics, they were assumed to be located in a sheltered environment and assigned the minimum wind−wave exposure value. Several environmental variables showed some collinearity (Pearson correlation r = 0.61−0.71 between mud, TN, and TP); however, this collinearity was within limits acceptable for tree-based machine learning methods such as GF (r < 0.9; Elith et al. 2010, Dormann et al. 2013).

Relative importance of environmental variables for predicting compositional turnover
GF was used to investigate estuarine benthic macro invertebrate turnover in response to land-derived stressors and natural environmental variables . Incidental taxa (≤ 3 occurrences across the entire dataset, n = 34) were not included in the GF models. The GF model had 2 components: the production of RF models (Breiman 2001) for each of the 88 input taxa using the R package 'extendedForest' (Liaw & Wiener 2002) and the aggregation of the individual split points from these models to calculate species turnover along each environmental gradient using the R package 'gradientForest' . Separate RF models were run for each species. Each RF model describes the relationship be tween an individual taxon and environmental variables by fitting an ensemble of regression models (1000 in our study) that each recursively split the observations into partitions. The proportion of out-of-bag variance explained measures the predictive power of the individual RF models (R 2 f ) , and the importance of each environmental variable (R 2 ) is measured by quantifying the degradation in performance when each environmental variable was randomly permuted ). This R 2 value described by Pitcher et al. (2012) and Ellis et al. (2012) refers to a unitless measure of cumulative importance. It should not be confused with the more commonly used R-squared (R 2 ) denoting coefficient of determination.
The degree to which abundance changes across adjoining partitions in each RF model represents a measure of compositional turnover occurring at the split value . GF aggregates the values of the tree splits from the RF models for all taxa models with positive fits (R 2 f > 0) to construct nonlinear empirical functions that represent predicted compositional change along each environmental gradient for the entire assemblage , hereafter referred to as compositional turnover. The compositional turnover function is measured in dimensionless R 2 units, where species with highly predictive random forest models (high R 2 f values) have a greater influence on the turnover functions than those with low predictive power (lower R 2 f ). The shapes of these monotonic turnover curves describe the predicted rate of compositional change along each environmental predictor; steep parts of the curve indicate fast assemblage turnover, and flatter parts of the curve indicate more homogeneous regions .
We extended the GF approach by adding a measure of uncertainty to the compositional turnover functions by bootstrapping GF models 100 times, similar to other regression tree-based methods (Leathwick et al. 2006). That is, the macroinvertebrate dataset was randomly sampled (with replacement) for each bootstrap iteration. The bootstrapping process was repeated 100 times, and at each iteration, compositional turnover functions were used to transform the environmental layers. Mean ± SD estimates of taxa R 2 f and environmental variable importance (R 2 ) were calculated for each GF model from the 100 bootstrapped iterations. To examine which taxa characterized the compositional turnover along each environmental gradient, cumulative abundance changes for the 5 taxa that achieved the highest cumulative importance values across the entire environmental gradient were also plotted.
Compositional turnover for each environmental predictor was visualized using principal component analysis (using the function 'prcomp' in the R package 'stats') to provide a multidimensional representation of variation in inferred community composition. Environmental variables were overlaid as vectors, indicating the strength and direction of the most important variables. All statistical analyses were undertaken in the software R (version 3.4.3; R Core Team 2019).

Relationships between compositional turnover and macroinvertebrate diversity
GLMs were used to explore the relative importance of compositional turnover along land-derived stressor gradients (mud, TN, TP) and natural environmental gradients (SST, SOI, wind−wave exposure) in explaining patterns in species richness (S; the number of taxa), Pielou's evenness (J' : Pielou 1966), Shannon-Wiener diversity (H';Shannon 1948), and numbers of rare taxa (those occurring only once or twice for each site/time). These 4 parameters will be referred to collectively as diversity indices. To be consistent with the GF models, incidental taxa (≤ 3 occurrences across the entire dataset, n = 34) were only included when calculating rare taxa (not in calculations of S, J', and H'). For the GLMs, the outputs of the GF model (compositional turnover values along 6 environmental gradients) were used as predictor variables, with compositional turnover along natural environmental gradients accounting for spatial and temporal dependency in the models. These compositional turnover values can also be thought of as environmental predictors that have been transformed to better reflect the underlying biodiversity patterns. This is because the turnover functions transform the measurement units of each environmental predictor to common biological units of compositional turnover. Data exploration was carried out according to the protocol developed by Zuur et al. (2010). Collinearity among predictor variables was generally low (Pearson's r < 0.5), with moderate correlations found only between turnover associated with TN and TP (r = 0.55), mud and TP (r = 0.66), and mud and TN (r = 0.74). The lack of strong correlations, and variance inflation factor values < 3, indicated that regressive models (including GF) should be able to separate land-based and natural variation (Zuur et al. 2010).
Models were fitted with error structures appropriate for the distribution of the data using the 'stats' and 'glmmTMB' (Brooks et al. 2017) packages in the software R (version 3.6.1; R Core Team 2019). A Poisson distribution with a log link function was used to model S and the number of rare taxa, a beta distribution with a logit link function was used for J', and a Gaussian distribution with an identity link function was used for H'. Interactions between predictors were already accounted for in the GF analysis; therefore, no interactions were included in the GLMs. Parsimonious models were developed using backward selection based on Akaike's information criterion (AIC) values to determine which variables were important in predicting patterns in diversity indices. As compositional turnover values were on the same scale, the relative importance of landderived stressors and natural environmental variables in predicting patterns in diversity indices was assessed using regression coefficients, with standard errors used as a measure of uncertainty. Model assumptions were verified by plotting residuals against fitted values, against each covariate in the model and against geographical coordinates (Zuur & Ieno 2016). Final models were checked for stability by varying the order in which variables were removed.

Relative importance of environmental variables for predicting compositional turnover
On average, across the 100 bootstrapped model runs, GF was able to effectively model species turnover for 82 ± 0.02 (SD) of the 88 input taxa (mean R 2 f = 0.49 ± 0.04). Both natural environmental variables and land-derived stressors were important in predicting patterns of compositional turnover in estuarine benthic macroinvertebrate communities, with the 3 natural variables combined slightly more important (27% of the conditional importance) than the 3 landderived stressors combined (22% of the conditional importance) overall (Fig. 2). SST (mean R 2 = 0.10) and wind−wave exposure (mean R 2 = 0.10) had the greatest influence on compositional turnover, followed by TP (mean R 2 = 0.08) and mud (mean R 2 = 0.08), TN (mean R 2 = 0.07), and SOI (mean R 2 = 0.06).
Non-linear patterns in compositional turnover were observed across all environmental gradients,  (Fig. 3). Sections of rapid turnover (steep sections of the curve) were observed along the wind−wave exposure, TP, mud, and TN gradients, indicative of large changes in species abundance and composition, followed by a levelling off indicating more homogeneous communities. For SST, high rates were initially followed by a slowing until 20°C and a rapid increase thereafter. The variability in mean predicted cumulative changes in compositional turnover, measured by the 95% prediction intervals, was relatively low. This uncertainty differed between environmental predictors and was greatest for TP and TN and lowest for SOI. Uncertainty also varied along individual predictor gradients with greater uncertainty observed where fewer data were available to inform predictions (higher levels of mud, TN, TP, and exposure gradients). Taxa identified as being important in characterising compositional turnover differed between the environmental gradients, although some taxa were characteristic of 2 or 3 environmental variables Fig. 3 Nicon aestuariensis along the exposure gradient). Rapid changes in abundance were generally associated with low variability, measured by the 95% prediction intervals, while higher variability was associated with flatter parts of the curves. Note that di rectionality of taxa response cannot be determined from these plots as they represent cumulative changes in abundance, that is, changes could be either increases or decreases in abundance at a given point along the gradient. Using these compositional turnover functions, shifts in community composition along environmental gradients were visualized in multivariate space where coordinate position represents inferred biological community composition, as associated with the environmental predictor variables (Fig. 5). The first 2 axes of the ordination plot captured 68% of the total variance. They demonstrated that both natural environmental variables (SST, wind−wave exposure, and SOI) and land-derived stressors (mud, TN, TP) were important variables influencing compositional turnover. Land-derived stressors influenced compo-sitional turnover in a similar way and, along with wind−wave exposure and SOI, were important in explaining biodiversity patterns along the first PC axis. SST also had a strong influence on compositional turnover, with site/times along the second PC axis showing high correspondence with location along the north to south gradient of New Zealand.

Relationships between compositional turnover and macroinvertebrate diversity
GLMs were used to determine whether compositional turnover driven by land-derived stressors and natural environmental variables resulted in positive or negative effects on macroinvertebrate diversity. The GLMs explained 7.8−13.4% of the variation in diversity indices, and all of the variables retained in the models were significant (p < 0.05), except for TP (p = 0.052) and wind−wave exposure (p = 0.121) in the model for H' (Table A1 in the Appendix). As hypothesized, compositional turnover along land-derived stressor gradients was linked to lower species richness, evenness, and diversity, and fewer rare taxa (Fig. 6). Compositional turnover along the sediment TN gradient had a negative effect on all 4 diversity , wind−wave exposure) derived from a gradient forest model. Points closer together indicate similarities in inferred community composition between site/times, and colour provides an indication of the latitudinal gradient (see Fig. 1 indices and was greater than the effect of turnover along other land-derived stressor gradients. Compositional turnover associated with increasing sediment mud content was only important in explaining patterns in S, while turnover associated with increasing sediment TP was only important in explaining patterns in H'. Compositional turnover along the SST and wind−wave exposure gradients had a positive effect on predicted values of J', H', and the number of rare taxa, with SST having a slightly stronger effect. Turnover along the SOI gradient was not important in explaining predicted patterns for any of the diversity indices. Greater uncertainty in model predictions was associated with compositional turnover along the TN and TP gradients (coefficient SE 2.0−3.5) compared with turnover along mud (coefficient SE 1.2) or natural environmental gradients (coefficient SE 0.8−1.5; Fig. 6, Table A1).

DISCUSSION
We have demonstrated that both land-derived stressors and natural environmental variables were important predictors of compositional turnover in estuarine benthic macroinvertebrate communities across New Zealand, reflecting a matrix of processes operating across multiple spatio-temporal scales. As expected, compositional turnover along land-derived stressor gradients was negatively associated with macroinvertebrate diversity indices, while turnover along natural environmental gradients (increasing SST and wind−wave exposure) generally had a positive relationship with these values.

Compositional turnover along natural environmental gradients
Predictably, across our large study area with its complex ocean currents influenced by both warm tropical and cold Antarctic water (Carter 2001), SST was the most important variable influencing compositional turnover. Temperature is known to be a key factor structuring communities across broad geographic scales (Tittensor et al. 2010, Denis-Roy et al. 2020, despite natural habitat characteristics such as grain size and salinity being important on local scales (Engle & Summers 1999, Denis-Roy et al. 2020. In our study, high rates of compositional turnover occurred above 20°C SST, corresponding to samples from the far north and east of New Zealand, where ocean temperatures have been increasing over the past 3 to 4 decades (Schiel 2013, Sutton & Bowen 2019. This high turnover rate suggests that climate change could lead to large shifts in community composition as physiological temperature tolerances are reached and species distributional boundaries change (e.g. Southward et al. 1995, Sagarin et al. 1999, Johnson et al. 2011. It is unlikely that temperature is the only driver of this compositional turnover pattern, however, with potential for it to act as a surrogate for a range of unmeasured broad-scale variables operating across the latitudinal gradient (e.g. species dispersal patterns, water circulation patterns, seasonality; Hawkins 2001, Thrush et al. 2005). Indeed, latitude was found to be an important driver of spatial patterns in fish assemblages across New Zealand (Stephenson et al. 2018), and a general latitudinal gradient in beta diversity has been observed in global-scale studies (Hillebrand 2004, Soininen et al. 2007, Qian et al. 2009), with higher species turnover toward the equator. These latitudinal patterns may arise because the physical limiting factors or ecological and evolutionary processes that influence turnover are also affected by latitude (Qian et al. 2009). In our study, compositional turnover along the SST gradient had a positive relationship with J', H', and the number of rare taxa, but not S. The pattern suggests that compositional turnover alters the relative proportion of rare to common species along this gradient, with common species becoming rarer with increasing temperature. Thus, the number of rare species increases with turnover associated with increasing SST, but the total number of taxa does not.
Wind−wave exposure, another important driver of species distributions in estuarine environments (Warwick et al. 1991, Hewitt & Thrush 2009, Hewitt et al. 2016, Denis-Roy et al. 2020, was the next most important variable influencing compositional turnover in this study. Although exposure and mud content often co-vary, these variables were not highly correlated in this study (r = −0.2), and the GF model would have accounted for any interactions between these 2 variables, suggesting no confounding effect. This is further supported by the different taxa characterising turnover along the mud and exposure gradients, with the only shared species (the crab Austrohelice crassa) exhibiting dissimilar changes in abundance along the 2 gradients. Like SST, compositional turnover along the wind−wave exposure gradient had a positive relationship with J', H', and numbers of rare taxa. In this study, sites with high exposure were located on central sandflats of a particular estuary type (large shallow drowned valley estuaries), for which the fetch allows wind-generated circulation and mixing. These high-energy environments generally have lower rates of sediment deposition, greater potential for recovery following storm events (Norkko et al. 2002b, Thrush et al. 2003a, improved food supply (via increased organic seston flux and/or resuspension of particulate organic matter; Fréchette & Bourget 1985, de Jonge & van den Bergs 1987, and increased potential for recruitment (Commito et al. 1995), which may promote diverse benthic communities in these areas.
Of the 6 environmental variables considered in this study, SOI explained the least amount of variation in compositional turnover, and turnover along the SOI gradient was not important in explaining patterns in estuarine benthic diversity. SOI influences a range of potentially important drivers (e.g. wind, temperatures, water column productivity) that could affect population dynamics and has been shown to be an important predictor of the abundance of species and functional traits (Hewitt & Thrush 2009, Hewitt et al. 2016. Unlike the other environmental factors considered in this study, SOI is a large-scale phenomenon that predominantly varies in time rather than space. The lack of robust time-series data for many sites in our analysis may have reduced the importance of this variable in predicting patterns of turnover compared with spatially variable factors.

Compositional turnover along land-derived stressor gradients
In our study, land-derived stressors were less important than SST and wind−wave exposure in predicting compositional turnover patterns in estuarine benthic communities across New Zealand. This result suggests that natural environmental variables regulate species distributions, with land-derived stressors constrained to act upon these existing communities. Given the low levels of mud and nutrients across many of our study sites, which are representative of estuaries across New Zealand, we would not expect land-derived stressors to be the most important variables influencing compositional turnover on a national scale. However, it is unknown whether the relative importance of the environmental variables would change if this model was applied to a dataset where levels of land-derived stressors were consistently high.
Once mud and nutrient levels were high enough to start acting as stressors to benthic communities, they began to have a discernible effect on compositional turnover despite the influence of natural environmental variables. High rates of turnover were observed between 0 and 10% mud, consistent with multiple studies that have shown a decline in functional redundancy and the abundance of sensitive taxa once mud content reaches 5−10% (e.g. Thrush et al. 2003b, Anderson 2008, Robertson et al. 2015, Ellis et al. 2017. For example, taxa characterising turnover along the mud gradient included the clam Austrovenus stutchburyi and the polychaetes Aonides spp., which showed rapid changes in abundance between 0 and 10% mud; these species have known preferences for sandy sediments with less than 10% mud (e.g. Norkko et al. 2002a, Gibbs & Hewitt 2004, Anderson 2008, Ellis et al. 2017). In contrast, the more constant changes in the abundance of the polychaetes Scolecolepides spp. and the mud crab Austrohelice crassa may reflect the tolerance of these species for a wider range of sediment grain sizes (e.g. Thrush et al. 2003b, Ellis et al. 2006, Anderson 2008, Robertson et al. 2015. In our study, high rates of compositional turnover in response to nutrients were observed at 1200 mg kg −1 , similar to a threshold of 1000 mg kg −1 TN associated with a shallowing of apparent redox potential discontinuity (aRPD) to near zero depths in 8 Californian estuaries (Sutula et al. 2014). Shallowing of the aRPD is usually associated with hypoxic events, which can lead to reduced abundance and diversity of benthic macroinvertebrates. Rapid changes in the abundance of specific taxa were observed at lower levels of nutrients, demonstrating that management thresholds based on compositional turnover will not protect all species. For example, rapid changes in the abundance of the bivalves Zemysia spp. and the polychaetes Magelona spp. were observed between 100 and 400 mg kg −1 TN, which is reasonably consistent with predicted distributions of these species between 200 and 600 and between 300 and 550 mg kg −1 TN, respectively (Ellis et al. 2017). The plateauing of compositional turnover observed at 3250 mg kg −1 TN and 1100 mg kg −1 TP, nutrient levels indicative of polluted estuaries (e.g. Oviatt et al. 1984, Gillespie & MacKenzie 1990, Sánchez-Moyano et al. 2010, may reflect a loss of taxa as communities become dominated by a limited number of species tolerant of high enrichment. Indeed, the GLMs showed that compositional turnover along these nutrient gradients was associated with lower species richness and diversity. However, the wide prediction intervals associated with these TN and TP thresholds mean these values should be interpreted with caution, as fewer data were available to inform the model. These values are reported as a contribution to the literature on nutrient effects and should be used in a weight-of-evidence approach in combination with other information, rather than relied upon as strict thresholds of community change along enrichment gradients. Additional sampling targeting locations with high levels of nutrients, as well as comparisons with thresholds identified using other ap proaches (e.g. Threshold Indicator Taxa ANalysis, Ecosystem Interaction Networks; Baker & King 2010, Thrush et al. 2021), would build confidence in the generality of these critical transitions.
Consistent with our second hypothesis, compositional turnover along land-derived stressor gradients was generally associated with lower S, J', and H' values, and fewer rare taxa. Maintaining diversity is important for promoting stability and resistance to disturbance (Levin et al. 2001), while rare taxa can confer functional resilience and make disproportionately large contributions to community and ecosystem functioning (Ellingsen et al. 2007). The loss of rare species has been proposed as an early warning signal of ecological shifts and functional impairment associated with anthropogenic stress as more of the community becomes represented by fewer, tolerant taxa . Across our study sites, compositional turnover along the sediment TN gradient was more important in explaining patterns in diversity indices than turnover associated with mud or TP, possibly because nitrogen is often the limiting nutrient in coastal systems (Howarth & Marino 2006). For example, compositional turnover along the nitrogen gradient could be linked to eutrophication-driven species loss. However, the importance of compositional turnover driven by both TN and TP in explaining patterns in H' suggests these nutrients can affect diversity in different ways. Similarly, turnover along the sediment mud gradient was also important in explaining patterns in S, highlighting the influence that multiple stressors can have on benthic diversity. The distinct groups of taxa characterising each of the land-derived stressor gradients also supports the idea that these stressors affect community turnover in different ways.
Hydrodynamic controls on sedimentation rates and nutrient loading can naturally result in upper reaches of estuaries being muddier and more enriched than outer reaches. While we cannot definitively conclude that human activities were the cause of elevated mud and nutrient levels in our study, we have shown that compositional turnover along these environmental gradients results in benthic macroinvertebrate communities with lower species richness, evenness, and diversity and fewer rare taxa. We also observed rapid changes in the abundance of functionally important species, such as Austrovenus stutchburyi and Macomona liliana, along land-derived stressor gradients. These bivalves influence community structure and microphytobenthic productivity as well as a range of physical and biogeochemical processes (e.g. sediment stability, pore water oxygen concentrations, nutrient cycling; Lelieveld et al. 2004, Thrush et al. 2006, Sandwell et al. 2009, Volkenborn et al. 2012. Consequently, these are the changes likely to occur if the total area of an estuary classified as being muddy or nutrient-enriched expands, with notable follow-on effects to ecosystem functioning (e.g. macroinvertebrate-mediated nutrient cycling; Lohrer et al. 2010) and, ultimately, the ecosystem services upon which humans rely. With increasing pressure on land worldwide, these land-derived stressors are likely to become more persistent. Even without intensification of human impact, the frequency and intensity of rainfall and storms are predicted to increase with climate change, likely increasing sedimentation rates and nutrient loading in estuaries (Inman & Jenkins 1999, McLean et al. 2001).

Consideration of uncertainty
Failure to consider uncertainty can result in poor management decisions (Regan et al. 2005, Link et al. 2012. Accordingly, we extended the GF approach by adding a measure of uncertainty to the compositional turnover functions and the changes in the cumulative abundance of key taxa. This development allowed results to be presented as an average of what is likely, thereby reducing the influence of non-representative outcomes. Indeed, for a single model run, SOI was the third most important variable explaining variation in compositional turnover, but averaged across 100 model runs, its relative importance decreased. Consistent with Sultana et al. (2020), who found that evenness of the environmental gradient can affect GF model performance, variability estimates associated with the compositional turnover functions and the changes in the cumulative abundance of key taxa indicated greater uncertainty where fewer data were available to inform predictions. For the key taxa, however, high rates of change were associated with low variability, providing confidence in these estimates.
Uncertainty also varied between environmental variables, with slightly less confidence associated with predictions of compositional turnover along TN and TP gradients. Although not explored explicitly in our study, greater uncertainty associated with nutrients could indicate that compositional turnover in response to nutrient loading is context-dependent. For example, high turnover may occur when nutrient loading and warm temperatures coincide, fuelling primary production, but a different response may occur if the same level of nutrient loading takes place in winter. Uncertainty may also be influenced by the restricted distribution of key taxa characterising these gradients, which may reflect habitat preference or sampling bias.
The addition of uncertainty estimates into GF outputs has important implications for management, which are not fully explored in this paper. For example, results could be spatially mapped (e.g. , Stephenson et al. 2018, Couce et al. 2020, with accompanying maps of uncertainty, to show the distribution of benthic communities and the uncertainty associated with those predictions. In this study, we might expect maps to highlight greater levels of uncertainty related to predictions of communities influenced by high levels of nutrient loading. Uncertainty was considered in the GLMs by comparing the size of the standard errors. Like the GF model, there was greater uncertainty linked to compositional turnover values along the TN and TP gradients in terms of predicting patterns in diversity indices.

Conclusions
We have demonstrated that both land-derived stressors and natural environmental variables, operating across multiple spatio-temporal scales, shape patterns of compositional turnover in estuarine macroinvertebrate communities across New Zealand. In this study, GF enabled us to tease out the effects of land-derived stressors from natural variation and identify critical levels where compositional turnover was high. Using GLMs, these turnover values were linked to measures of benthic macroinvertebrate diversity, which indicated that turnover along landderived stressor gradients had a negative effect on benthic communities at a local scale. Relationships identified by these exploratory models are correlative, and while they do not necessarily prove a causal link, they do identify possible drivers of patterns that could be investigated further through controlled experiments . Exploratory models also allow for studies to be undertaken on much larger scales than funding for manipulated experiments would allow, providing information about pro-cesses operating over broad scales. Future work could examine other environmental variables, including biotic factors (e.g. competition for resources, predation, small-scale biological disturbance), incorporate measures of environmental variability (e.g. seasonal ranges of predictors rather than averages), and consider lag effects. GF also allows for the inclusion of abundance data from different survey methods ) because a dimensionless R 2 measure is used to quantity compositional turnover, meaning that other estuarine taxa, such as fish, could be included in models to provide a more holistic view of ecosystem response. This study moves towards an ecosystem-based management approach by considering how multiple land-derived stressors influence patterns of estuarine compositional turnover and diversity, against a background of natural variability operating at multiple spatio-temporal scales.
Acknowledgements. We thank Massey University and the New Zealand Ministry of Business, Innovation and Employment (MBIE) for funding this work through the Oranga Taiao Oranga Tangata research programme (contract MAUX1502, led by Murray Patterson). Additional funding was provided by the Cawthron Institute's Internal Investment Fund and the National Institute of Water and Atmospheric Research Coasts and Oceans Programme (Project: COME2001). We also acknowledge the support of New Zealand regional authorities that provided data and permission to use it: Bay of Plenty Regional Council, Environment Canterbury, Environment Southland, Greater Wellington Regional Council, Hawkes Bay Regional Council, Marlborough District Council, Nelson City Council, Northland Regional Council, Otago Regional Council, Tasman District Council, and West Coast Regional Council. Paula Casanovas (Cawthron Institute) provided helpful advice throughout the writing of this manuscript. We also thank 4 anonymous reviewers for constructive comments that improved the manuscript.