Site characteristics

This study was undertaken in the Manawatu River catchment, and all sites were on the mainstem or tributaries of the Manawatu River (Fig. 1). One or 2 sites at each of 7 rivers (Makakahi, Manawatu, Mangatainoka, Oroua, Oruakeretaki, Tiraumea, and Tokomaru) were selected for Phormidium monitoring (Table 1, Fig. 1). Sites were chosen to represent rivers with and without historical Phormidium proliferations, with similar geology, and a range of nutrient and flow regimes (Table 1). Land cover was primarily heavy pastoral or native vegetation, with urban accounting for <1% in catchments (Table 1). Riparian or overhead shading from vegetation was negligible at all sites. Sites 2/3 and 6/7 were paired sites within 1.9 and 0.62 km of each other but situated down- and upstream (respectively) of sewage treatment plant (STP) discharges (Fig. 1). Surveying and sampling were undertaken approximately weekly at 7 sites (sites 1, 4–7, 9, and 10) from 6 January 2012 to 26 or 28 June 2013, at 2 sites (sites 2 and 3) from 18 July 2012 to 28 June 2013, and at site 8 from 17 August 2012 to 28 June 2013.

Site surveys and sampling

Surveys were conducted in runs and riffles. The length of the reach surveyed varied from 30 to 100 m among sites. At each site, 4 transects perpendicular to the water’s edge and extending to a maximum depth of 0.6 m were surveyed. Percent Phormidium cover was estimated visually at 5 points along each transect with the aid of an underwater viewer (giving a total of 20 views/site; Ministry for the Environment and Ministry of Health 2009). Substrate composition (macrophytes, moss, mud [<0.2 mm], sand [0.2–2 mm], fine gravel [2–8 mm], gravel [8–64 mm], cobble [64–264 mm], boulder [>264 mm], and bedrock) was assessed visually at the same 20 viewpoints.

Water temperature, conductivity, and pH were measured once at each site on each sampling occasion using a handheld multiparameter PCTESTR 35 sonde (Oakton Instruments, Vernon Hills, Illinois), and water samples were collected for dissolved nutrient analysis. When Phormidium mats were present, a single sample was collected from 10 rocks at each site by scraping mat material from a defined circular area (6.25-cm diameter). The 10 samples were pooled, chilled, and stored in darkness for transport to the laboratory where they were stored frozen (−20°C; within 5 h) until toxin analysis.

Continuous river flow data were obtained from permanent flow-gaging stations in close proximity to each sampling site (Fig. 1). Daily mean flow was calculated for each site. Long-term median flows at each site were obtained from Horizons Regional Council (Table 1).

Laboratory analysis

Nutrients in filtered water were analyzed on a Lachat QuickChem^® Flow Injection Analyser (FIA+ 8000 Series, Zellweger Analytics, Loveland, Colorado) based on the methods provided by APHA (2005). The nutrients analyzed were: NH₄-N, NO₂-N, NO₃-N, and dissolved reactive P (DRP). The accredited detection limits were 0.01 mg/L NH₄-N, 0.002 mg/L NO₂-N and NO₃-N, and 0.005 mg/L DRP.

Phormidium mat samples were lyophilized and homogenized. Lyophilized material (100 mg) was suspended in 10 mL of Milli-Q water (MQ; Millipore, Billerica, Massachusetts) containing 0.1% formic acid and sonicated (30 min; Cole Parmer 8890; Biolab, Auckland, New Zealand) on ice. Samples were centrifuged (3000g, 10 min). The supernatants were transferred to a septum-capped vial and analyzed for ATX, HTX, dhATX and dhHTX based on liquid chromatography-mass spectrometry (LC-MS) as described by Wood et al. (2016). External standards were used for ATX calibration, using dilutions of a certified reference material for ATX (National Research Council, Canada) to prepare working standards that ranged from 0.5 to 20 ng/mL in 0.1% formic acid. A relative response factor of 1, with ATX as the calibration reference, was used to quantify HTX, dhATX, and dhHTX.

Data analyses

Generalized additive mixed models (GAMMs; Hastie and Tibshirani 1990) were used to model nonlinear trends in % Phormidium cover and total toxin concentration (sum of ATX, HTX, dhATX, and dhHTX in mg/kg of lyophilized material) in relation to time (wk) of the year and environmental variables. The % Phormidium cover data were means from 20 views obtained at each approximately weekly sampling occasion (n = 677). Total toxin models were fitted only with observations from sampling occasions when Phormidium samples were collected for toxin analysis (n = 176). Both % Phormidium cover and total toxin data were strictly positive, highly right-skewed, over-dispersed, and with a large proportion of 0 values, so models were constructed with log-normal errors. Seven noncollinear predictor variables were selected on the basis of initial data exploration and ecological knowledge and were included in the models. Collinearity among predictor variables was checked using the variance inflation factor (VIF < 3; Zuur et al. 2010). Week of the year was included as a continuous covariate (1–52) to account for seasonal trends based on a cyclic cubic spline to allow the continuity between weeks 1 and 52. Biggs and Close (1989) suggested that simple point-by-point correlations between nutrient concentrations and periphyton biomass do not provide a true indication of the water-column nutrients to which periphyton has been exposed over its accrual period. To address this problem, cumulative mean DIN (the sum of NH₄-N, NO₂-N, and NO₃-N) and DRP concentrations over the accrual period were calculated for each time point (i.e., the mean concentrations since the previous sampling point at which the % Phormidium cover was 0). ‘Accrual’ DIN and DRP were √(x)- and log(x)-transformed, respectively. The influence of flow was considered by including log(x)-transformed ‘times median flow’, i.e., the maximum flow magnitude (obtained from daily mean flow data) between sampling periods divided by the long-term median flow for a given site. Water temperature, accrual DIN, accrual DRP, times median flow, and conductivity were included as continuous covariates, and site as a categorical variable with 10 levels. Percent Phormidium cover also was included as a continuous covariate in the total toxin model. To account for temporal autocorrelation in the weekly time-series data, models were fitted based on autoregressive moving average correlation structure (ARMA) of order 1. Models were selected with a stepwise procedure based on the generalized Akaike Information Criteria (GAIC) and were validated by inspecting the deviance residuals. Final models are presented as partial effects plots, which show the effect of each predictor variable conditional to others in the model. The partial effects of each predictor are displayed as cubic splines showing either negative or positive effects relative to the overall mean of the response variable centered on 0. Partial plots also show standard errors around the fitted spline and partial residuals for each observation. The nonlinear interaction between accrual DIN and DRP on % Phormidium cover was further investigated with a tensor product smoother in a GAMM model. The interactive effect of DIN and DRP on % Phormidium cover was visualized in a contour plot of the predicted values.

We also investigated the magnitude of river flow required to reduce % Phormidium cover at each site to <20%. Analysis was undertaken only for the 4 sites that had % Phormidium cover >20% on multiple occasions (sites 1, 3, 5, and 10). For each site, the maximum flow (obtained from daily means) divided by the long-term median flow on the 7 d prior to surveying was selected. A time frame of 7 d was selected because this is the length of time over which % Phormidium cover has been observed to increase from 0% to >20% cover (Heath et al. 2011, Wood and Young 2012). Greater than 20% Phormidium cover was selected because this level is the lower threshold of the amber alert level in the New Zealand Guidelines for Cyanobacteria in Recreational Freshwater (Ministry for the Environment and Ministry of Health 2009). If a site has a cover of 20%, cover probably will continue to increase in the absence of a flushing event. For each site, quantile regression models (80^th, 85^th, and 90^th percentiles) were fitted to the % Phormidium cover data with log(x)-transformed river flow as the predictor. The river flow at which % Phormidium cover stayed at 0% was estimated from the regression lines.

All statistical analyses were performed with the software R (version 3.1.1; R Project for Statistical Computing, Vienna, Austria), GAMMs models were run with the gamlss package (Rigby and Stasinopoulos 2005), and quantile regressions were run with the quantreg package (Koenker 2012).

Substrate composition, spatial and temporal patterns in % Phormidium cover, physicochemical variables

The dominant substrate at all sampling sites was large (12–25 cm) and small (6–12 cm) cobbles (>42% at all sites; Fig. S1). Phormidium mats were present at all 10 river sites, but mean cover was spatially and temporally variable (Fig. 2). Five sites (1, 2, 5, 6, and 10) had mean Phormidium cover >40% on multiple occasions. Three sites (7, 8, and 9) had high % cover in single views (40–60%), but the mean cover was ≤15% (Fig. 2). At 2 sites (3 and 4) mean cover was <2% (Fig. 2).

At the sites with high % Phormidium cover, the period with the greatest % cover was generally from early austral summer (November) to the end of autumn (May). Periods of high cover varied among sites over the study period. For example, high cover was recorded at site 5 and 10 from mid-April to mid-June 2013, but this pattern was not observed at other sites (Fig. 2). Although situated <2 km apart, % Phormidium cover differed notably between sites 2 and 3 and sites 6 and 7. In general, sites downstream of STPs had higher cover. This situation was particularly apparent at site 2 in March 2013 and site 6 in January 2013 (Fig. 2).

Median DIN concentrations varied markedly among sites. The lowest concentration was recorded at site 4 (0.12 mg/L) and the highest at site 6 (0.87 mg/L; Fig. 3A). Median DRP concentrations were generally low at all sites (0.01 mg/L). The 3 highest concentrations were recorded at sites 2, 6, and 7 (Fig. 3B). Median water temperatures ranged between 9.5 and 15.5°C and were comparable across sites except site 4, which was cooler (Fig. 3C). The lowest (4.6°C) and highest (26.2°C) water temperatures were recorded at site 4 in July 2012 and site 5 in February 2013, respectively. Median conductivity was variable among sites (range: 53–313 µS/cm at sites 4 and 9, respectively; Fig. 3D).

Percent Phormidium cover in relation to environmental variables

The most parsimonious GAMM explained 56.2% of the total deviance in % Phormidium cover and included all predictor variables (Table 2; p < 0.05). Week of the year affected predicted % Phormidium cover positively for the first half of the year (January–May; p < 0.001; Fig. 4A). The model predicted greater cover with increasing accrual DIN concentrations up to 0.8 mg/L, after which the effect remained constant (p < 0.001; Fig. 4B). Accrual DRP had a weak, marginally significant negative effect on % Phormidium cover when DRP > 0.005 mg/L (p = 0.03; Fig. 4C). Water temperatures >15°C had a strongly positive effect, whereas colder temperatures did not affect % Phormidium cover (p < 0.001; Table 2, Fig. 4D). Times median flow had a strongly positive effect for flows <0.4 and a negative effect for flows >0.5 (p < 0.001; Fig. 4E). Percent Phormidium cover increased monotonically with conductivity (p < 0.001; Fig. 4F). Percent Phormidium cover differed strongly among sites (p < 0.05; Fig. 4G).

Sites could be grouped loosely into 3 categories based on the interactive effects of accrual DRP and DIN on % Phormidium cover (Fig. 5A): 1) low DRP (<0.01 mg/L) and elevated DIN (>0.1 mg/L; sites 1 and 5), 2) downstream of STPs with elevated DIN and DRP (>0.1 mg/L and 0.01 mg/L; sites 2 and 6), and 3) low DRP (<0.01 mg/L) and DIN (<0.06 mg/L; site 10). The GAMMs predicted % Phormidium cover >20% when DIN was between 0.09 and 0.27 mg/L and DRP was between 0.004 and 0.005 mg/L (Fig. 5B).

Predicting the flow required to reduce % Phormidium cover to <20%

Quantile regression based on the 90^th percentile of the maximum daily mean flow (divided by the long-term median) between the sampling periods predicted that the flow required to reduce % Phormidium cover to <20% would vary markedly among the sites 1, 3, 5, and 10 (~2.5, ~0.7, ~5.4, and ~249× times the long-term median, respectively; Fig. 6A–D). At site 10, the magnitude of required flushing flow predicted based on the 90^th percentile was large, so we also explored the 80^th and 85^th percentiles. These regressions predicted required flushing flows that were ~6.75 and ~16.5× the long-term median, respectively, at site 10 (Fig. 6D).

Anatoxin concentrations and variant composition

Anatoxins were detected in variable concentrations in samples from 6 of the 10 sites (Fig. 7A, B, Table S1). The highest median toxin concentration was measured at site 10 (0.20 mg/kg). The highest toxin concentration measured in a single sample was at site 1 (82.4 mg/kg, 17 February 2012). When % Phormidium cover was >~10%, the relationship between cover and toxin concentrations was close to linear (Fig. 7B).

No ATX variant was detected in any sample (Fig. 7C). Over the entire sampling period, HTX was the dominant variant at 4 sites (2, 8, 9, and 10), where it accounted for >50% of the total toxin concentration (Fig. 7C). dhHTX was dominant at 2 sites (1 and 5; Fig. 7C). The dominant variant showed some temporal variability among sites, but no obvious patterns were observed (Table S1).

Total anatoxin concentration and relationship with physicochemical variables and % Phormidium cover

The most parsimonious GAMM for total toxins explained 50.8% of the total deviance and included 5 predictor variables: week of the year, accrual DIN, times median flow, conductivity, and % Phormidium cover (Fig. 8A–E). The effect of week of the year was weak and marginally significant (p = 0.06; Fig. 8A), and followed a similar pattern to that observed for the % Phormidium cover GAMMs. The effects of accrual DIN and flow were less pronounced than those for % Phormidium cover, but were significant (p < 0.001), with the highest toxin concentrations predicted at mid-range DIN and during low flows (Fig. 8B, C). Low conductivity ~50 µS/cm) resulted in high toxin concentrations. Higher conductivity produced a similar pattern but had a smaller effect on toxin concentrations (Fig. 8D). Percent Phormidium cover had a strongly positive and almost linear effect on toxin concentration (p < 0.001; Fig. 8E).

Drivers of % Phormidium proliferations

The P. autumnale-complex (Strunecký et al. 2010) is one of the most widespread cyanobacterial taxa in the world and has been identified in a variety of habitats and nutrient regimes including eutrophic and oligotrophic ponds, lakes, and rivers, and on damp soils (Komárek and Anagnostidis 2005, Richter et al. 2009, Wood et al. 2016). Nevertheless, it and other Phormidium species have received international scientific attention only recently because of their ability to form toxin-containing benthic proliferations (Gugger et al. 2005, Bouma-Gregson and Higgins 2015, McAllister et al. 2016). Research on environmental drivers of Phormidium proliferations in cobble-bed rivers has been spatially and temporally limited but has highlighted the importance of variables including river flow, fine sediment, and water-column nutrients (Heath et al. 2011, 2015, Hart et al. 2013, Wood et al. 2015, McAllister et al. 2016). Our study was limited in its geographic extent, but spanned the range of river water-chemistry conditions observed in New Zealand rivers that experience Phormidium proliferations. The GAMMs indicated that week of the year, water chemistry, temperature, and flow are important regulators of Phormidium proliferations.

The basic model for control of periphyton biomass in cobble-bed rivers identifies hydrologic disturbance as the primary regulator, within which nutrients operate by influencing the rate of accrual during periods of stable flow (Biggs 1995, Biggs et al. 2005). Our findings differ from those of other investigators who showed that Phormidium was present in higher abundance than other algal species when DRP concentrations were elevated (Horner et al. 1990), was prevalent in rivers with high DIN and DRP concentrations (Loza et al. 2013), and was out-competed by other cyanobacteria when grown in competitive assays under low-nutrient regimes (Loza et al. 2014). However, these investigators did not identify Phormidium to species level, and different species and strains within genera may vary in their optimal nutrient requirements (Heath et al. 2016). Moreover, the mats studied by Horner et al. (1990) and Loza et al. (2013, 2014) contained many cyanobacterial or algal species. In mats similar to those we observed, Phormidium was volumetrically abundant (Brasell et al. 2015) and formed thick and cohesive mats. This growth structure probably is critical to understanding the prevalence of Phormidium mats in rivers with low-nutrient concentrations.

Sand-Jensen (1983) noted that the relationships between water-column nutrients and algal mats are most important and most evident during early growth stages before mats mature and internal processes, such as nutrient recycling, begin. Several mechanisms might underlie the ability of Phormidium to form expansive mats under low-nutrient conditions; e.g., formation of geochemical conditions (e.g., pH and dissolved O₂) within mats that enable release of DRP bound to sediments entrapped in the mat (Wood et al. 2015), internal microbial processes (i.e., N fixation by bacteria in the mats), and possible advection of DIN and DRP from turbulent near-bed flow that could compensate for low-nutrient concentrations. Future investigators should explore nutrient concentrations in Phormidium mats directly to provide insights to this possible disconnection between nutrients in the water column and actual nutrient availability.

Previous investigators have suggested that water-column DIN concentrations >0.10 to 0.20 mg/L promote Phormidium proliferations (Wood and Young 2012, Heath et al. 2015, McAllister et al. 2016). This threshold held true for 4 of the 5 sites in our study that experience Phormidium proliferation. However, DIN concentrations at site 10 are much lower than the 0.10 mg/L threshold, but the site experiences prolonged proliferations (Fig. 2). Reasons for this phenomenon are unknown but may be related to the higher intensity of flow required to remove Phormidium mats at this site than at other sites (see Discussion below), which would provide longer periods for mats to become established or would result in larger starting inoculums after flushing events. These relic patches might already contain functioning microbial communities, e.g., N-fixing organisms, that would enable rapid mat expansion.

Conductivity had a strongly positive effect on % Phormidium cover. Conductivity is often correlated with periphyton biomass and is strongly related to catchment geology and water source in New Zealand rivers (Biggs 1990, 1995). It also is linked integrally to river flow, and generally increases during periods of low flow because of lack of dilution from water with lower ionic content, a trend also evident in our data set (data not shown). Many of the major ions contributing to conductivity, e.g., Ca²⁺, Cl⁻, Na⁺, K⁺, and Mg²⁺, influence important metabolic processes or enhance/suppress growth in cyanobacteria (Seale et al. 1987, Parker et al. 1997, Carneiro et al. 2011). Individual ions/compounds that contribute to conductivity may play a role in regulating Phormidium accrual, but further research is required to test this suggestion.

Our analysis highlighted the importance of week of the year, with higher % Phormidium cover predicted in summer and early autumn than at other times. At most sites, this period coincides with periods of stable flow and elevated temperature, variables that also were related to elevated cover. Our model predicted greater % Phormidium cover when water temperature was >15°C. A recent analysis of a 20-y data set of benthic algae in 2 pristine Norwegian streams also forecast greater Phormidium cover with increasing temperatures (Schneider 2015). River flow equal to approximately half the long-term median flow was projected to favor % Phormidium cover, with a ‘hump-back curve’ apparent (Fig. 4E). A habitat suitability criteria model developed for Phormidium based on velocity, depth, and substrate type showed a similar nearly parabolic curve, with highest cover observed at moderate velocities (1 m/s; Heath et al. 2015). Gas bubbles caused by photosynthesis are often visible in Phormidium mats. Bubble formation is more likely to occur under low flow because diffusion of O₂ will be slowed by the thick boundary layer (Hawes et al. 2014). Therefore, Phormidium growth in slower flow may be coupled with more frequent autogenic detachment, which would reduce biomass accrual in low flows (Boulêtreau et al. 2006). Thus, flows that are sufficient to enhance nutrient and gas flux but insufficient to generate shear stress that results in biomass loss are postulated to be optimal for Phormidium accrual (sensu Biggs and Thomsen 1995).

Intensity of flushing flow required to reduce Phormidium proliferations

High-energy flow events ‘reset’ periphyton communities (Clausen and Biggs 1997). Some management authorities involved in issuing human health warnings related to Phormidium have used ‘length of time since a flushing flow’ as an early warning indicator of elevated risk. For example, Milne and Watts (2007, p. 24) advocated for the criteria “…that if there has not been a flushing flow (defined in their study as 3× the median flow [sic]) for two weeks and low river flows (set at the lowest 10^th percentile flow for each river).” Exceeding these criteria triggered weekly surveys of sites prone to Phormidium proliferations.

We used quantile regression to predict the magnitude of flushing flow required to reduce % Phormidium cover and showed that this value differed among sites. Identifying the reasons for these differences is challenging. Phormidium mats are likely to be removed during periods of elevated flow through high shear stresses, abrasion by mobilized sediments, and grinding action of tumbling gravel/cobble substrata (Grimm and Fisher 1989, Horner et al. 1990, Biggs et al. 1999, Francoeur and Biggs 2006). Heath et al. (2015) showed in an intensive field survey that the best predictor for Phormidium presence was stable substrates unlikely to be disturbed by flood flows. Substrate heterogeneity also might be relevant if cobbles or boulders provide more refuges (i.e. cracks and crevices) during flushing flows and enable faster recolonization (Bergey 2005, Murdock and Dodds 2007). In addition, Phormidium mats undergo a series of developmental stages as they grow. These stages include initial attachment, maturation, and dispersal (McAllister et al. 2016). Dispersal can occur when bubbles of O₂ become trapped in the mats and cause them to detach or slough off the substrate. Early-stage mats are more firmly adhered to substrates and, therefore, are likely to require more intense flushing to remove them.

Anatoxins

Total anatoxin concentrations and the relative composition of structural variants showed high spatial and temporal variability among and within sites, an observation consistent with previous studies (Heath et al. 2011, Wood et al. 2010). Toxin concentrations were generally low compared to those in a nationwide data set (McAllister et al. 2016), and the complete absence of toxins in mats from sites 3, 6, and 7 is noteworthy.

An important caveat when assessing relationships between toxin concentrations in Phormidium mats and physicochemical variables is that mats contain a mixture of toxic and nontoxic genotypes (Heath et al. 2010, Wood et al. 2010), other organisms (Brasell et al. 2015), and inorganic debris (Wood et al. 2015). These variables can vary among sites. Moreover, even within toxic genotypes, the relative amount of toxin produced can vary 100-fold (Wood et al. 2012). This variability limits our ability to infer relationships among physicochemical variables and actual toxin production or shifts in genotypes. However, our analysis provides some indication of periods or conditions associated with greatest toxin content of mat samples and, therefore, periods of highest risk to river users.

We identified 3 prominent relationships in our analysis. First, total toxin concentrations were higher when conductivity was low than when it was high, a pattern opposite that observed for % Phormidium cover. A similar result was noted in a culture-based study where the anatoxin quota was lowest in high N and P conditions (Heath et al. 2016). The authors suggested that this result supports the growth–differentiation balance hypothesis, which states that actively dividing cells are less likely to produce secondary metabolites than are nondividing cells. Second, the relationship between % Phormidium cover and toxin concentrations was nearly linear, particularly when % Phormidium cover was >10% (Figs 7A–C, 8A–E). This relationship differs from findings of other studies that suggested no relationship between these variables (Wood et al. 2010, Heath et al. 2011) and culture-based studies that showed peak toxin quotas during the initial growth phase (Harland et al. 2013, Heath et al. 2016). Third, toxin concentrations were lower when DIN was <0.2 mg/L, a pattern similar to that observed for % cover and consistent with the observations by Heath et al. (2016) but not with the growth–differentiation balance hypothesis. Development of molecular techniques (e.g., quantitative polymerase chain reaction) that enable enumeration of Phormidium cells in environmental samples and distinguishing toxic and nontoxic genotypes is essential to allow insights into how these and other variables are related to toxin production and genotype succession. This ability may ultimately assist in understanding the ecological role of anatoxins in these mat communities.