Catchment land use predicts benthic vegetation in small estuaries

Statistical analysis

We first screened predictors for high collinearity. If there were sets of predictors with pair-wise correlations >0.7, we eliminated all predictors bar one. The retained predictor was the one with the lowest sum of correlations with predictors other than those in the inter-correlated set. Land use and riverine DIN concentrations were highly correlated (r = 0.84), so DIN concentration was excluded from the analysis because it is dependent upon land use. We scrutinized the distributions of the retained predictors. Several were extremely right-skewed, so these were log-transformed (designated by † in Table 1). Once the distributions were near normally or near uniformly distributed, we standardized (mean = 0, standard deviation = 1) the predictors to make the ranges of all predictors comparable and to assist in model convergence.

We used two approaches to identify the potentially important predictors and to identify the relative importance of these predictor variables. First, we used Bayesian variable selection using stochastic search (O’Hara & Sillanpää, 2009). This method identifies those predictors that had high posterior probabilities of being included in the best models for explaining variation in the MA:TV ratio. We used the posterior odds ratio framework to assess predictor importance (Kass & Raftery, 1995). A predictor is assigned an uninformative prior for being included in the best model (a predictor is equally likely to be selected as not), which corresponds to a prior odds ratio of 0.5 (included):0.5 (not included) = 1. If the posterior probability of inclusion, after calculations, is (or exceeds) 0.75, then the posterior odds are 0.75 (included):0.25 (not included) = 3. The ratio of the posterior odds to the prior odds is the posterior odds ratio (here 0.75:0.25/0.5:0.5 = 3), with values exceeding 3 being indicative of probable importance of a predictor in explaining variation in the response variable (here MA:TV ratio). Models were calculated using JAGS (Plummer, 2003).

Second, we used hierarchical partitioning (HP) on the predictor variables to calculate the relative proportions independently explained by each predictor. We used the hier.part (Walsh & Mac Nally, 2004) package in R (R Development Core Team, 2011). HP complements Bayesian model selection by quantifying the relative amounts of variation independently explained by each predictor (Mac Nally, 1996).

The %fertilized predictor proved to be important (Table 1) but we were concerned that its effect might be moderated by the catchment to estuary (C:E) ratio or residence time of the estuary (Tf). Therefore, we used Bayesian model selection and HP analyses for a full interaction model involving these three predictors, notwithstanding that the C:E ratio and Tf were not important from the analyses in Table 1. The full interaction model included the three predictors, each pair of interactions, and the three-way interaction (see Table 2).

We fitted a change-point model for the relationship between MA:TV and %catchment fertilization (F). The model was: ${\displaystyle MA:T{V}_i\sim Normal\left({\mu }_i,\sigma \right),}$ ${\displaystyle {\mu }_i=\alpha \ast \delta \left(F_i-\gamma \right)+{\beta }_1\ast \delta \left(\gamma -F_i\right)\ast F_i+{\beta }_2\ast \delta \left(F_i-\gamma \right)\ast F_i;}$

where: δ is unity if the argument is non-negative and zero otherwise and γ is the change-point. The priors were: $\alpha ,{\beta }_i\sim Normal\left(0,\sigma =2\right),\mathrm{and}\gamma \sim Uniform\left(0,100\right)$. The program JAGS (Plummer, 2003), which uses a Gibbs sampler, was used to fit the relationship; there were 12,000 iterations and 5,000 ‘burns-in’ samples. We checked convergence using Gelman–Rubin methods (convergence of multiple independent chains).

Land use as a predictor of shallow benthic vegetation

The results show that land use is a strong predictor of the proportion of macroalgae to total vegetation within south-eastern Australian estuaries. Although the current analysis shows that land use is a stronger predictor than nitrogen loads, we do not interpret this to mean that nitrogen inputs to estuaries are not the key driver of changes to estuarine benthic vegetation. Rather we use these findings to shed light on the possible mechanisms through which nutrients drive change within estuaries and how catchment land use integrates this change.

The proportion of the catchment receiving fertilization was a better predictor of the MA:TV ratio than was areal dissolved inorganic nitrogen (DIN) load, which we have previously suggested to be a better predictor of MA:TV in these estuaries than total nitrogen (TN) or total phosphorous (TP) loads (Woodland et al., 2015). This outcome arose because there was a relatively weak relationship between measured loads (both total and normalized to the estuary area) and the total area fertilized (km²) within catchments (R² = 0.33). The lack of such a relationship is consistent with previous studies that have shown nitrogen attenuation factors can be highly variable (Elwan et al., 2015). Therefore, the relationship between land use and estuarine response is not just driven by a land-use-load relationship, as we had expected.

There was a strong non-linear relationship between DIN concentration in the rivers and the MA:TV ratio (R² = 0.78, when DIN concentration is log transformed) arising from the relationship between land use and DIN concentrations within the rivers (Woodland et al., 2015). However, we do not believe that the nutrient concentrations observed within the rivers are the primary driver of changes in the MA:TV ratio because these rivers drain into estuaries of greatly different sizes and hence dilution. Moreover, there was no relationship between DIN and the MA:TV ratio when the nutrient concentration was normalized to estuary area. These results suggest that catchment land-use metrics ‘integrate’ factors affecting the amount and availability of nutrients within the estuary that control the MA:TV ratio, which are missed by instantaneous measurements of load. Catchment land-use metrics may incorporate: (1) the historical sequence of delivery of nitrogen (N) and total suspended solids that are trapped and recycled or re-suspended within the estuary; (2) increased bioavailability of particulate and dissolved organic N delivered to estuaries as fertilization increases (Seitzinger, Sanders & Styles, 2002; Petrone, Richards & Grierson, 2009); and (3) local groundwater inputs of N directly to estuaries (Wong et al., 2014).

Our results suggest that estuarine vegetation structure can be substantially altered when agricultural land use constitutes as little as 24% of the catchment. Therefore, it is instructive to compare the nutrient loading and export rates measured here with previous studies to place the land-use intensity in this study in a wider context. The areal loading rates of N in these estuaries span the range reported globally for estuaries, ranging 10²–10⁵ mmol N m⁻² yr⁻¹ as total N (Woodland et al., 2015). The rates of N generation from the catchments in the lowest quartile of %fertilization averaged 0.9 and 0.25 kg ha⁻¹ yr⁻¹ for TN and DIN respectively (Table 3), which are at the lower end of DIN exports of 0 to 10 kg ha⁻¹ yr⁻¹ reported in forested catchments (Bernal, Butturini & Sabater, 2005; Brookshire et al., 2012). Our undisturbed catchments have lower exports than forested catchments elsewhere in the world, highlighting the relatively oligotrophic state of estuaries fed by pristine catchments in Australia. For the most fertilized catchments (>80% fertilization), we saw average N generation rates of ∼4.5 and 1.9 kg ha⁻¹ yr⁻¹ for TN and DIN respectively, which are comparable with reported nutrient generation rates for mixed farming/rural land use in southeastern Australia (Drewry et al., 2006). Our N generation rates are at the lower end of reported nutrient generation rates of 4–14 kg ha⁻¹ yr⁻¹ for TN in European and North American systems (Howarth et al., 1996), highlighting that even small amounts of relatively low-intensity agriculture can lead to large changes in benthic vegetation in these naturally oligotrophic estuaries. Studies from other locations are needed to investigate whether the patterns observed here are globally applicable.

The use of the macroalgae to total vegetation ratio as an indicator of estuarine condition

It is virtually impossible to select an ecological indicator that represents all critical aspects of ecosystem function. Our choice of MA:TV was based on the requirement that we could easily obtain relevant data for large areas of the estuary. In shallow estuaries, such as those studied here, macroalgae is widely considered an indicator of eutrophication (Valiela et al., 1997). There are cascading ecological consequences from the increasing dominance of macroalgal biomass to food webs, from changes in consumer biodiversity, productivity and trophic relationships (e.g., omnivory) to biogeochemical cycling and dissolved-oxygen dynamics (Sogard & Able, 1991; Valiela et al., 1997). Consistent with this, we also saw that once catchment fertilization exceeded 24%, and alongside a transition to macroalgal dominance of demersal vegetation, all chlorophyll a measurements were >6 µg/L, which typically is regarded as eutrophic (Hakanson, Bryhn & Hytteborn, 2007).

The ability to reliably predict MA:TV ratio using just one variable differs from previous studies that have shown that multiple predictors are needed to explain >50% of the variation of other response variables (Li et al., 2007; Greene et al., 2015). One of the strongest relationships from previous reports has been between chlorophyll a and the area of agriculture land use in a catchment and estuary volume (R² of 0.68) in 14 Canadian estuaries (Meeuwig, 1999). This is not unexpected because a certain nutrient load will be diluted to different extents depending on estuary volume. Similarly, one would expect phytoplankton concentration to be sensitive to estuarine residence time, which will lead to different wash-out rates (Nixon et al., 2001). We saw no clear relationship between chlorophyll a and percent fertilization in our data set, which was consistent with the need for other variables to describe the response of this parameter.

We found no important relationship between seagrass or total vegetation areal extent and land use. Elsewhere, a combination of land-use and physical factors, such as tidal range and mean wave height, were needed to describe seagrass areal extent (Li et al., 2007), illustrating the interplay of factors other than eutrophication in controlling seagrass distribution. By standardizing macroalgal extent as a proportion of total vegetation, our analysis reduces the influence of physical factors, such as sediment movement and hydrodynamics, that often limit the growth of benthic vegetation. The MA:TV measure also accounts for estuaries with different hypsometric profiles because the MA:TV ratio is functionally constrained to those areas where light penetration can support benthic vegetation. Change point analysis showed that the MA:TV ratio increased at a lower rate above a catchment fertilization of 24%. This probably suggests that any increase in biomass above this point may have manifested itself as increased thickness (as opposed to area). Alternatively, macroalgae became growth limited at high biomass due to limitation by other factors such as light and/or space. Therefore, a disadvantage of this approach is that it does not sensitively distinguish between moderate and high levels of disturbance.

The interaction model showed that estuary flushing time did not contribute much to explaining variation in the MA:TV response. The residence times used in our study are relatively short (0.6–4.2 days compared to months to years for lagoons), which may partially explain the lack of importance of residence time. However, macroalgae can assimilate N in several hours and, given the subsequent relatively slow turnover of N, residence time may not significantly affect the macroalgal response (Nixon et al., 2001). Lagoon systems, with much longer residence times, are likely to respond differently to our estuaries because phytoplankton are more dominant in systems where water residence time exceeds phytoplankton turnover time (Hauxwell & Valiela, 2004). As the catchment-estuary area ratio (C:E ratio) increased, we expected that nutrient inputs would be distributed over a smaller estuarine area and may render estuaries more sensitive to the proportion of fertilizing land uses in the catchment. The results of the interaction model, which included C:E ratio (Table 2), suggested that there was little evidence of an interaction of the C:E ratio with the proportion of fertilizing land uses in the catchment. As the C:E ratio increased, the transit time of loads delivered to the estuary decreased, leading to lower retention and exposure to nutrients within the system compared to estuaries with lower C:E ratios. Systems with low C:E ratios have the load spread over a larger area, but with a longer transit time, leading to higher retention and exposure to nutrients. Our results suggest that these opposing effects may largely cancel each other out, leading to the C:E ratio having no perceptible effect on MA:TV.

Application to management

Although nutrient loads are a critical management tool for receiving waters, the expense and long time frame required to collect and meaningfully interpret these data mean that such data are not always available. However, land-use information typically is much more readily available, and as has been illustrated here and previously (Meeuwig, 1999; Meeuwig, Kauppila & Pitkänen, 2000), can provide a good indicator of likely risk to estuaries. Once estuaries at risk have been identified, there should be a further assessment of ecological impact. The MA:TV ratio provides a relatively rapid and spatially representative indication of ecological response and condition. By way of an example case study, an Index of Estuary Condition (IEC) is being used to assist the prioritization of estuary management investment thereby supporting the Victorian Waterway Management Program, Australia (DEPI, 2013). Indicators that have demonstrated relationships with processes threatening estuaries (e.g., nutrient loading and land-use change) are essential if broad scale resource condition assessments are to be interpretable, ecologically meaningful and useful for management (Barbour et al., 2000; Stoddard et al., 2008). For these reasons the MA:TV index is being incorporated into the Victorian IEC.