Impact of wastewater treatment upgrade and nitrogen removal on bacterial communities and their interactions in eutrophic prairie streams

Abstract

Eutrophication can impact bacteria by altering fluxes and processing of nutrients and organic matter. However, relatively little is known of how bacterial communities, diversity, and interactions with phytoplankton might respond to nutrient management. We used 16S rRNA amplicon sequencing to compare bacterial assemblages in the water column upstream (control) and downstream (impact) of a wastewater treatment plant (WWTP) located on a eutrophic prairie stream. Sampling occurred before (2012) and after (2018) the 2016 biological nutrient removal (BNR) upgrade that removed >90% of nitrogen (N, mainly NH₄⁺). Multivariate ordination suggested that effluent-impacted bacterial communities were associated mainly with elevated NH₄⁺ concentrations before the upgrade, whereas those after BNR were characteristic of reference systems (low NO₃⁻, diverse regulation). Genera such as Betaproteobacteria and Rhodocyclacea were abundant at impacted sites in 2012, whereas Flavobacterium and a potential pathogen (Legionella) were common at impacted sites in 2018. Nitrifier bacteria (Nitrospira and Nitrosomonas) were present but rare at all sites in 2012, but recorded only downstream of the WWTP in 2018. Generalized additive models showed that BNR reduced bacterial diversity, with ∼70% of the deviance in diversity explained by hydrology, pH, nutrients, and phytoplankton abundance. Overall, NH₄⁺ removal reduced symptoms of cultural eutrophication in microbe assemblages.

16S rRNA, ammonia, biological nutrient removal, effluent, phytoplankton, Qu’Appelle River

Introduction

Urban and agricultural development has intensified eutrophication and negatively impacted water quality worldwide through increases in nitrogen (N) and phosphorus (P) fluxes to surface waters (Carpenter et al. 1998, Preisner et al. 2020). Enrichment of aquatic ecosystems causes excessive production of algae, harmful cyanobacterial blooms, and deposition of organic matter, which leads to oxygen depletion, fish kills, and alterations to the biogeochemistry of metals and other sparingly soluble compounds (Dodds and Smith 2016; Jane et al. 2021). Under natural conditions, productive streams and rivers have high abundances of diverse microbes, with particularly high levels of heterotrophic bacteria relative to algae and cyanobacteria (Bernhardt et al. 2018, Hamdhani et al. 2020). Such lotic ecosystems act as ‘active pipes’ that process organic matter and nutrients along the hydrological flow path (Peterson et al. 2001). In contrast, anthropogenic eutrophication of flowing waters can disrupt multiple ecosystem functions (biotic interactions, organic material processing, and nutrient cycling) and human services (recreation, fishing, and potable water) while converting flowing waters into less reactive conduits (Peterson et al. 2001, Dodds et al. 2009, Hamdhani et al. 2020). Under these conditions, factors regulating the interaction between heterotrophic and autotrophic assemblages are less well understood (Bernhardt et al. 2018, Hamdhani et al. 2020), and there remains a need to understand how bacterial communities respond to changes in nutrient pollution and management.

Advanced wastewater treatment processes are designed to prevent the degradation of aquatic ecosystems and re-establish natural functions in receiving freshwaters (Tchobanoglous et al. 2003). Treatment of urban septic waters often involves sequential primary, secondary, and tertiary processing, which removes particulate organic matter (POM), dissolved organic matter (DOM), and dissolved phosphorus (P), respectively, but not total dissolved nitrogen (TDN) or ammonium, NH₄⁺ (Hamdhani et al. 2020). As a result, tertiary wastewater processing is often a substantial point source of ammonium to downstream habitats (Hamdhani et al. 2020). More recently, biological nutrient removal (BNR) has used sequential nitrification and denitrification by microbes to remove N, especially NH₄⁺, while also providing effective DOM and P control (Carey and Migliaccio 2009). Consequently, application of BNR technologies or other N-removal processes (e.g. ion exchange, reverse osmosis, adsorption, advanced oxidation process, chemical precipitation, and air stripping) to urban wastewater treatment plants (WWTP) can result in sudden declines in N pollution (Hamdhani et al. 2020, Kabuba et al. 2022) and provide an opportunity to study the environmental impacts of selective N removal on bacterial communities (Bergbusch et al. 2021a,b). However, to date, fewer than 30% of cities employ BNR to remove N from urban effluent, even in developed countries (Tchobanoglous et al. 2003, Hamdhani et al. 2020), in part because of scientific and management concerns about the detrimental effects of selective N removal and low N:P ratios on eutrophication of downstream freshwaters (Schindler et al. 2016).

Heterotrophic bacterial communities in freshwater ecosystems exhibit variation in diversity and relative abundance in response to changes in nutrient pollution. For example, the diversity of benthic bacterial communities initially declines in response to the influx of nutrients and antimicrobial compounds (Drury et al. 2013), although diversity can recover downstream if new taxa are introduced by effluent or if the degree of pollution is limited (Wakelin et al. 2008, Price et al. 2018). In general, the influx of NH₄⁺-rich effluent from tertiary WWTPs promotes Sphingobacteriales and ammonia-oxidizing bacteria such as Nitrospira and Nitrosomonas (Drury et al. 2013, Sonthiphand et al. 2013). In contrast, the influx of BNR-treated effluent is known to favour growth of orders of the phylum Bacteroidetes, Proteobacteria, and Chloroflexi that prefer oxic, less nutrient-rich environments (Atashgahi et al. 2015, Price et al. 2018). Together, these patterns suggest that the transition from tertiary treatment (P and DOM removed) to BNR (N, P, and DOM removed) may reveal the unique effects of NH₄⁺ on the structure and diversity of bacterial communities in recipient ecosystems, although, as yet, no study has tested this hypothesis.

Changes in urban effluent flux and composition may also alter relationships between heterotrophic bacteria and photoautotrophs (as biofilms or phytoplankton) in receiving freshwaters (Bernhardt et al. 2018). For example, the growth of select bacterial phyla (Proteobacteria, Bacteroidetes, Actinobacteria, Firmicutes, and Deinococcus-Thermus) in rivers, lakes, and the ocean is correlated negatively or positively with the extant abundance of cyanobacteria, indicating a mutualistic relationship based on nutrient syntrophy (Berg et al. 2009). Similarly, inputs of nutrients and pharmaceuticals from urban effluent (Oliveira and Goulder 2006, Waiser et al. 2011) can limit the nutrient exchange between planktonic bacteria, algae, and cyanobacteria, thereby altering community composition and potentially efficacy of aquatic DOM processing (Sonthiphand et al. 2013, Price et al. 2018). Finally, while selective removal of NH₄⁺ from urban wastewater can reduce primary production, favour algae over cyanobacteria, and affect plankton-benthic interactions (Bergbusch et al. 2021a,b), little is known of the consequences of these changes to the structure and function of bacterial consortia, especially in lotic ecosystems. Improved understanding of how bacteria, algae, and cyanobacteria interact in response to changes in nutrient status is essential for the development of effective surface water management strategies (Carpenter et al. 1998).

In this paper, we employed 16S rRNA amplicon sequencing to assess bacterial communities in two lotic ecosystems and understand the effects of tertiary and BNR-treated effluent on planktonic bacterial assemblages in an N-limited eutrophic stream continuum. First, we quantified the effects of nutrient pollution on the composition of planktonic bacteria, algae, and cyanobacteria in water samples collected at nine stations above and below of WWTP outfall from the City of Regina, Saskatchewan, Canada, both before (2012) and after (2018) the 2016 transition from tertiary to BNR wastewater treatment. Second, we used generalized additive models (GAMs) and redundancy analysis (RDA) to evaluate the relative importance of biological interactions and physical-chemical conditions as the mechanisms controlling the relative abundance and diversity of suspended bacterial communities. Third, we sought to isolate the unique effects on bacterial assemblages of a >90% decline in pollution with NH₄⁺, but not P or DOM (∼0% change) due to WWTP upgrade to BNR processes. We hypothesized that: (a) bacterial community composition at wastewater-impacted sites would be different from upstream or control sites in both years; (b) there would be a shift from eutrophic (nutrient-adapted and cycling taxa) to more oligotrophic taxa after BNR; and; (c) bacterial diversity would decrease due to a reduction in N and phytoplankton after commissioning of BNR processes (Atashgahi et al. 2015, Price et al. 2018). We also predicted that bacterial composition and diversity would be affected by declines in cyanobacterial and chlorophytes in favour of diatoms (Bergbusch et al. 2021b), separate from any effects of variation in the physical-chemical environment (Kent et al. 2007).

Methods

Sampling area

The Qu’Appelle River drainage basin is located in southern Saskatchewan, Canada (50.77°N, 103.99°W), mainly within Treaty Four territory (Fig. 1). The area has a subhumid continental climate with annual precipitation of ∼380 mm, most rainfall occurring between May and July, and runoff from snowmelt during the spring (Vogt et al. 2018). Qu’Appelle River and Wascana Creek are P-rich, turbid, low-order (1–3) grassland streams (Leavitt et al. 2006, Bergbusch et al. 2021a,b). Within the basin, the Qu’Appelle River arises near Eyebrow Lake and flows eastward through Buffalo Pound Lake, Pasqua Lake, and five other lakes before draining into the Assiniboine River. Wascana Creek originates south of the City of Regina from groundwater and flows northwest before joining the Qu’Appelle River mid-reach between Buffalo Pound and Pasqua lakes. South-flowing sub-saline Last Mountain Creek also drains into the Qu’Appelle River upstream of Pasqua Lake, unless subject to hydrological management (Fig. 1). More than 90% of the Wascana Creek (Swarbrick et al. 2022) and the Qu’Appelle River catchments (Hall et al. 1999, Hosseini et al. 2018) are used for dryland grain farming or livestock with diffuse nutrients inputs, while urban nutrients are contributed from the cities Moose Jaw and Regina (Leavitt et al. 2006, Bergbusch et al. 2021a,b). Prior to 2016, the City of Regina removed P and DOM via tertiary chemical treatment (Leavitt et al. 2006), whereas N (mainly NH₄⁺) was additionally removed thereafter through bacterial denitrification and anammox processes (Bergbusch et al. 2021a,b). The City of Moose Jaw uses a modified tertiary nutrient removal system in which some N was eliminated via denitrification in storage lagoons prior to diversion to agricultural fields or release into the Qu’Appelle River near Buffalo Pound Lake (Moose Jaw River Watershed Advisory Committees 2006).

Figure 1.

Map of Wascana Creek and the Qu’Appelle River in Saskatchewan, Canada. The study area is outlined in an inset of Canada.

We evaluated the effects of N removal by BNR on downstream bacterial communities in Wascana Creek and the Qu’Appelle River. Nine locations were sampled biweekly from May to September, both before (2012) and after (2018) the 2016 Regina WWTP upgrade to BNR (Fig. 1). The nine locations are distributed across 200 km of serpentine streams, with the first site being 50 km upstream of the WWTP, sites 3–5 being 3–60 km downstream, and the furthest site in the Qu’Appelle being 140 km from the wastewater plant. Two reference sites were located upstream of the WWTP in Wascana Creek, while three impacted locations were downstream of the WWTP prior to the confluence of Wascana Creek and the Qu’Appelle River. In the Qu’Appelle River, one site was located upstream of the confluence with Wascana Creek, whereas a further three sites were downstream of the confluence, prior to Pasqua Lake (Fig. 1).

Field sampling methods

Stations were sampled over two consecutive days every two weeks between May and September in 2012 and 2018 following Bergbusch et al. (2021a,b). Briefly, environmental parameters were measured using a YSI Model 85 meter for water temperature (°C), dissolved oxygen content (mg O₂ L⁻¹), pH, and specific conductivity (µS cm⁻¹) at surface and mid-column depths. Water discharge (m³ s⁻¹) was measured using a Swoffer Instruments Model 2100 current velocity meter. Water transparency (NTU) was recorded using a LaMotte 2020 model turbidity meter. Water samples were collected from each site at the river surface using distilled water-rinsed 10-L containers. These containers were simply submerged just below the surface at each site. This water was then screened through a 243-µm aperture mesh to remove coarse POM but not colonial cyanobacteria.

Water column laboratory analyses

Water from 2012 and 2018 was filtered onto a 1.2-µm glass-fibre filter (pigment analyses, frozen at −10ºC) or a 0.45-µm pore membrane filter (DNA analyses, frozen at −80ºC). Filtrate from the 0.45-µm pore filtration was frozen at ∼ −20ºC before analysis for nutrient chemistry. Dissolved ammonium (NH₄⁺/NH₃), nitrate (NO₂⁻ + NO₃⁻), and soluble reactive phosphorus (SRP; PO₄³⁻) were quantified (mg N L⁻¹ or µg P L⁻¹) using a Lachat QuikChem 8500 FIA automated ion analyzer following standard procedures (American Public Health Association, American Water Works Association and Water Environment Federation 1998).

Phytoplankton abundance and gross community composition were quantified by measuring changes in biomarker pigment concentrations using either trichromatic spectrophotometric analysis (chlorophyll a [Chl a], ug L⁻¹; total phototrophs) or high-pressure liquid chromatography (mainly carotenoids, nmoles L⁻¹; other algal and cyanobacterial taxa) following standard procedures (Leavitt and Hodgson 2001). Biomarker pigments were isolated and quantified on a model 1100 Agilent HPLC system fitted with a photodiode array and fluorescence detector. Detectors were calibrated with authentic standards following Bergbusch et al. (2021a,b). Biomarker pigments included total phototrophs (Chl a), siliceous algae and some dinoflagellates (fucoxanthin), mainly diatoms (diatoxanthin, diadinoxanthin), dinoflagellates (peridinin), chlorophytes (Chl b), cryptophytes (alloxanthin), colonial cyanobacteria (myxoxanthophyll), Nostocales cyanobacteria (canthaxanthin), potentially diazotrophic cyanobacteria (aphanizophyll), and total cyanobacteria (echinenone). However, most analyses were based on fucoxanthin, Chl b, and echinenone alone, as they represented the three predominant groups of suspended algae (‘phytoplankton’) in these streams (Bergbusch et al. 2021b).

16S rRNA gene sequencing and bioinformatic processing

POM on 0.45-µm pore-size membrane filters was frozen at −80°C and DNA extracted using a Qiagen Mo Bio Powersoil DNA isolation kit. Water-column samples from 2012 were initially frozen at −20°C until analysis, whereas samples from 2018 were filtered immediately and DNA extracted. DNA yield was quantified using the Qubit BR dsDNA assay (ThermoFisher Scientific) and diluted to 10 ng uL⁻¹. Libraries were prepared according to Kozich et al. (2013), using primers designed for the v4 hypervariable region of the bacterial S16S rRNA gene (515F, 5′-GTGCCAGCMGCCGCGGTAA-3′ and 806R, 5′-GG ACTACHVGGGTWTCTAAT-3′ from Wu et al. 2015). Sequencing of 16S rRNA genes was performed using an Illumina MiSeq Benchtop Sequencer, a MiSeq® Reagent Kit v2, and 500 cycles (2 × 250 bp paired-end).

Raw Illumina sequence reads were processed using the Mothur software (v.1.39.5; Kozich et al. 2013) to remove the low-quality and chimeric reads, eliminate artifacts, and reduce duplicity in the dataset (Schloss et al. 2009). DNA nucleotide sequences of the 16S rRNA gene were assigned to operational taxonomic units (OTU) through 97% identity alignment in the SILVA small subunit (SSU) reference bacterial database (Pruesse et al. 2007, version 132). All samples were subsampled to 1014 reads to allow for further statistical analysis.

Statistical analyses

All statistical analyses and visualization of sample data were performed using R version 4.0.4 software (R Core Team 2022). All data visualizations were created using ggplot2 v. 3.3.5 (Wickham 2016). Before modelling, we calculated Shannon Weiner diversity to evaluate bacterial community OTU richness and evenness responses to changes in physical and chemical conditions, as well as phytoplankton abundance and composition, along the stream continuum (Kim et al. 2017). Richness is the number of different taxa in a sample and evenness is the similarity of population sizes across taxa, and together these are calculated with the Shannon Weiner diversity formula: the negative sum of the proportion of the community represented by one OTU multiplied by the natural logarithm of this proportion (Kim et al. 2017). Other diversity metrics were considered, including species richness, Inverse Simpson, and Chao1 indices.

GAMs were used to quantify the linear and non-linear relationships between bacterial diversity, phytoplankton abundance and composition, and spatial and temporal variation in physical-chemical conditions (Pedersen et al. 2019). Models were generated using R packages mgcv v. 1. 8–39 (Wood 2011) and gratia v. 0.7.3 (Simpson 2022), with residual marginal likelihood smoothness selection (Wood 2011). Models were fitted with distance from the WWTP and calendar day of year (DOY) smooth terms to evaluate spatial and temporal variation in diversity, whereas other GAMs were developed with physical-chemical parameters (transformed for normalization as needed) as response variables to infer processes regulating bacterial diversity. Shannon diversity was normally distributed, so a Gaussian distribution with link identity function was employed, whereas discharge, Chl b, and echinenone were square-root transformed, and NH₄⁺, N, SRP, and fucoxanthin were log₁₀ transformed. All models included a ‘select = true’ feature to remove response variables that do not account for variation in the diversity data. Relationships between bacterial OTU diversity and predictor variables were based on a smooth and random effect developed by Wood (2013a,b). Ninety-five percent credible intervals (P = 0.05) were used to visualize uncertainty in significant physical-chemical variables. We assessed basis size, dispersion of residuals, homogeneity of variance, and the relationship between predicted and observed responses to ensure model assumptions were not violated.

Ordination models were created using the vegan v. 2.5.7 package in R (Oksanen et al. 2022) to summarize how the composition of suspended bacteria (at phylum level) varied across stream gradients (spatially), throughout time (seasonally), and in response to effluent influx. We used nonmetric multidimensional scaling (NMDS) to assign distances to sampled OTUs and visualize the similarity of bacterial communities. Distances were assigned using the Bray–Curtis dissimilarity index, with five dimensions selected to reduce stress to 0.1 and ensure goodness of fit. To complement the NMDS analysis, we did an analysis of similarities (ANOSIM) with a Bray–Curtis dissimilarity index to test for significant differences between site and year dissimilarity matrices.

RDA was selected to quantify how the composition of bacterial consortia varied in response to linear combinations of measured biological, chemical, and physical parameters. Discharge, Chl b (chlorophytes), and echinenone (cyanobacteria) were square-root transformed, whereas NH₄⁺, SRP, and fucoxanthin (diatoms) were log₁₀ transformed to normalize residuals. Predictors are represented as arrows within RDA plots, whereas samples are plotted as points within ordination space. Highly correlated predictors in RDA exhibit similar magnitude and direction, while communities that are correlated positively to the predictors plot near the arrow heads. Negative correlations plot parallel to the arrows, but in the opposite direction from the arrowhead.

Relative abundances of OTUs were plotted using the phyloseq library, Bioconductor version 3.12, and BiocManager 1.30.16) (Mcmurdie and Holmes 2013) in R to evaluate taxonomic changes that occurred following the BNR upgrade. We tested the enrichment of genera >2% of total abundance at sampling locations grouped together by treatment (control one, impacted sites, control two, and downstream sites) for each year with the LEFSe approach. The LEFSe is an approach used to identify genomic features (taxa, genes) that differ by condition (treatment type, location) based on a number of non-parametric tests, including the Kruskal–Wallis sum-rank test (determine effect of condition), Wilcoxon rank-sum test (determine biological consistency), and linear discriminant analysis (determine effect size) (Segata et al. 2011). To normalize genera abundance data to ensure representation by low values, we used the ‘copies per million’ option, which pre-sample normalizes the sum of the values to 1e+06.

Map creation

The map of the study area (Fig. 1) was created in QGIS (QGIS Development Team 2023) using layers from GeoGratis (Government of Canada 2023).

Results

Physical-chemical conditions

Nitrogen concentrations in Wascana Creek decreased up to 97% for NH₄⁺ and 33% for TDN following the WWTP upgrade to BNR, whereas there were few changes in P concentrations beyond those reflecting variation in background stream conditions (Table 1). Specifically, while BNR upgrade had relatively little effect on NO₃⁻ concentrations at sites downstream of wastewater outfall in Wascana Creek (‘Impacted’; Sites 3, 4, and 5) in both years, NH₄⁺ levels decreased sharply from 11.63 to 0.34 mg N L⁻¹ due to BNR. Similarly, TDP varied little among years downstream of urban effluent outfall, although SRP values decreased moderately from 0.34 to 0.07 mg P L⁻¹ reflecting background variation in SRP levels upstream of the WWTP (Site 2; Table 1). As a result, TDN:SRP ratios increased three-fold (30–95) in the wastewater-impacted sections of Wascana Creek and generally remained elevated in the downstream Qu’Appelle River after the BRN upgrade. Variation in the concentration of most nutrients was limited between years at control and post-confluence sites within the Qu’Appelle River (sites 6–9). Overall, we suggest that the strongest effects of effluent on the lotic system occurred in Wascana Creek within ∼60 km downstream of the wastewater outfall (Sites 3–5).

Table 1.

The mean and standard error of physical-chemical variables in Wascana Creek and Qu’Appelle River in 2012 and 2018.

Year	2012					2018
Location	Control 1	Control 2	Impact	Control 3	Downstream	Control 1	Control 2	Impact	Control 3	Downstream
	Wascana Creek			Qu’Appelle River		Wascana Creek			Qu’Appelle River
Stream site	1	2	3, 4, and 5	6	7, 8, and 9	1	2	3, 4, and 5	6	7, 8, and 9
TDN (mg N L−1)	1.19 ± 0.26	0.97 ± 0.32	6.85 ± 0.66	0.51 ± 0.06	1.47 ± 0.17	2.19 ± 0.18	1.09 ± 0.09	4.60 ± 0.44	0.69 ± 0.04	1.58 ± 0.19
TDP (mg P L−1)	0.15 ± 0.02	0.10 ± 0.02	0.16 ± 0.02	0.03 ± 0.01	0.05 ± 0.00	0.55 ± 0.09	0.09 ± 0.02	0.18 ± 0.03	0.02 ± 0.00	0.06 ± 0.01
SRP (mg P L−1)	0.28 ± 0.07	0.26 ± 0.04	0.34 ± 0.05	0.07 ± 0.01	0.12 ± 0.02	0.37 ± 0.07	0.05 ± 0.02	0.07 ± 0.01	0.01 ± 0.00	0.03 ± 0.01
TDN:SRP	8.53 ± 2.93	5.25 ± 2.67	30.40 ± 5.14	8.14 ± 0.83	16.21 ± 1.86	11.19 ± 4.25	48.46 ± 16.79	93.39 ± 12.78	61.07 ± 6.11	227.48 ± 146.59
Nitrate (mg N L−1)	0.34 ± 0.22	0.26 ± 0.12	1.65 ± 0.29	0.46 ± 0.26	1.78 ± 0.20	0.02 ± 0.01	0.04 ± 0.02	1.52 ± 0.19	0.02 ± 0.01	0.27 ± 0.07
Ammonium (mg N L−1)	0.07 ± 0.01	0.18 ± 0.06	11.63 ± 1.82	0.07 ± 0.03	0.27 ± 0.06	0.05 ± 0.01	0.09 ± 0.03	0.34 ± 0.12	0.03 ± 0.00	0.07 ± 0.02
Temperature (˚C)	17.64 ± 2.01	18.41 ± 1.87	18.80 ± 0.88	19.20 ± 1.47	19.49 ± 0.75	17.63 ± 1.25	17.41 ± 1.32	18.57 ± 0.6	17.73 ± 0.95	18.33 ± 0.56
DO (mg O2 L−1)	7.82 ± 0.71	6.70 ± 1.05	4.68 ± 0.63	6.40 ± 1.01	7.23 ± 0.23	7.28 ± 0.63	6.61 ± 0.76	10.03 ± 0.55	7.61 ± 0.38	9.75 ± 0.27
SPC (µS cm−1)	1840.38 ± 149.93	1104.63 ± 85.34	1569.67 ± 25.12	1095.63 ± 48.28	1414.57 ± 50.86	2212.40 ± 121.73	991.33 ± 62.87	1561.50 ± 25.87	872.31 ± 45.63	1409.93 ± 72.73
pH	9.04 ± 0.18	8.37 ± 0.10	7.78 ± 0.07	8.44 ± 0.09	8.45 ± 0.06	8.54 ± 0.21	8.60 ± 0.10	8.47 ± 0.09	8.30 ± 0.05	8.70 ± 0.07
Turbidity (NTU)	10.53 ± 2.96	25.30 ± 6.34	8.81 ± 1.84	62.71 ± 10.03	50.12 ± 4.93	6.53 ± 0.62	21.36 ± 2.29	10.66 ± 1.58	35.48 ± 5.35	25.42 ± 2.97
Discharge (m3 s−1)	0.79 ± 0.48	1.91 ± 1.26	2.26 ± 0.50	4.75 ± 1.68	10.16 ± 1.60	0.06 ± 0.05	0.53 ± 0.32	1.01 ± 0.13	1.62 ± 0.96	3.34 ± 0.38

TDN refers to total dissolved nitrogen, TDP refers to total dissolved phosphorus, SRP refers to soluble reactive phosphorus, and SPC refers to specific conductivity.

River discharge varied substantially among years at all sites, with two- to four-fold greater runoff in 2012 than in 2018 (Table 1). As a result, turbidity was slightly higher at reference sites in Wascana Creek and downstream in the Qu’Appelle River during the high flow of 2012. Overall, turbidity peaked immediately upstream of the WWTP (Site 2), then decreased at impacted sites for both sampling years before increasing again in the Qu’Appelle River. Despite variation in discharge, there were few marked differences among years in terms of temperature, specific conductance, pH, or dissolved O₂. However, dissolved O₂ was nearly two-fold greater after BNR than before the upgrade at effluent-impacted stations within Wascana Creek (Table 1).

Bacterial community composition

Ordination analysis with NMDS showed that bacterial communities were generally different between years at wastewater-impacted sites (Sites 3–5) but relatively similar at control (Sites 1, 2) and downstream (Sites 6–9) locations (Fig. 2A). Based on the analysis of similarity, both year (ANOSIM R = 0.258, P = 0.001) and site (ANOSIM R = 0.332, P = 0.001) factors were significant. Immediately downstream of wastewater input to Wascana Creek (Site 3), bacterial assemblages were similar between 2012 and 2018; however, communities at Sites 4 and 5 differed among years, with 2018 sites grouping outside of the overlap between the 95% confidence interval year ellipses (Fig. 2A). Prior to the WWTP upgrade, the bacterial assemblages at Sites 4 and 5 clustered near Site 3, suggesting a similarly impacted consortium at each station (Fig. 2B). In contrast, after BNR was instituted, communities at Sites 4 and 5 did not group with Site 3 and instead ordinated more tightly with the upstream Wascana Creek control site, suggesting some degree of recovery from urban nutrient pollution. Compositions of all bacterial assemblages in the Qu’Appelle River stations (Sites 6–9) during 2012 and 2018 appeared to be similar, clustering in the bottom left quadrant of the ordination (Fig. 2). Taken together, these findings suggest that differences in bacterial assemblage composition between 2012 and 2018 are generally more pronounced at wastewater impact sites (3–5) but not upstream or downstream sites.

Figure 2.

Non-metric multidimensional scaling for bacterial communities at each site in Wascana Creek (1–5) and Qu’Appelle River (6–9) in 2012 and 2018. Circles are 2012 communities and triangles are 2018 communities. Communities that group together are similar, whereas those located farther apart are more different. Distances were assigned by the Bray–Curtis dissimilarity index. The 95% confidence intervals are represented by ellipses, with a stress value of ∼0.1. In A, the ellipses are by year, whereas B ellipses are by site.

Changes in bacterial phyla and specific genera also reflected differences between sampling years, with some pronounced variation among sampling stations (Fig. 3, Supplementary Fig. S1). Most sites in both years had bacterial communities comprised of the phyla Actinobacteria, Bacteroidetes, and Proteobacteria with moderate increases in Proteobacteria and Bacteroidetes at impacted sites (3–5) in 2012 and 2018, respectively (Supplementary Fig. S1). At the genus level, some taxa exhibited enrichment during 2018 relative to 2012, although patterns were not necessarily consistent with the degree of wastewater pollution (Fig. 3). For example, Flavobacterium was notably more abundant at impacted sites, both relative to headwaters and in 2018 relative to 2012. In contrast, several taxa were more abundant in wastewater-impacted sites than in headwaters in either 2012 [Alcaligenaceae (unclassified), Alsobacter, Aquishaera, Betaproteobacteria (unclassified), Flavobacterium, Methyloparacoccus, Mycobacterium, Polynucleobacter, Rhizobiales (unclassified), Rhodocyclaceae (unclassified)] or 2018 [Betaproteobacteria (unclassified), Flavobaterium, Fluviicola, Hydrogenophaga, and Legionella]. This pattern suggests that specific variation in Flavobacterium spp. may be leading to observed increases in the phylum Bacteroidetes after BNR processes were commissioned (Figs 3 and  4), whereas Alcaligenaceae, Alsobacter, Betaproteobacteria, Polynucleobacter Rhizobiales, and Rhodocyclaceae were contributing mainly to changes in the phylum Proteobacteria.

Figure 3.

LEFSe bubble plots for water-column bacterial genera in control one (1–2), impacted (3–5), control two (6), and downstream (7–9) locations in Wascana Creek (1–5) and Qu’Appelle River (6–9) in 2012 and 2018. LEFSe is employed to determine if taxa at certain sites or treatments are significantly enriched relative to others. Only significant genera (P < 0.05) and those present at >2% relative abundance are displayed. Genera are associated with phyla by colour.

Figure 4.

Nitrifier genera abundance represented by the OTU 16S rRNA copy number at each site summed across sampling weeks.

Nitrifiers

Nitrifying bacteria were uncommon in water column samples in both years, although there were more detections across sites in 2012. Summed across sampling days from May to September, the number of 16S rRNA reads for nitrifiers (Nitrosomonadaceae [unclassified], Nitrosomonas, and Nitrospira) was low at between 0 and 3. Nitrifying taxa were detected at impacted sites in 2012 and 2018, while these bacteria were found only upstream and beyond heavily impacted sites in 2012 and at only one site (7) and one genus (Nitrospira) in 2018.

Response of bacteria community to physical-chemical conditions and phytoplankton

RDA suggested that bacterial communities exposed to wastewater (impact sites) were correlated strongly to variation in NH₄⁺ concentrations from tertiary-treated wastewater, whereas effluent-impacted assemblages were under more diverse control mechanisms after the BNR upgrade (Fig. 5). Specifically, bacterial communities at impacted sites (Sites 3–5) in 2012 ordinated along the second RDA axis in the left-most plot (20.1% variance explained) and were correlated positively with NH₄⁺ concentrations. In contrast, the impacted communities in 2018 loaded mainly on the first RDA axis (17.0% variation) in the right-most plot and were correlated positively with variation in NO₃⁻ and SRP. As well, to some degree, taxa at impacted Sites 4 and 5 grouped near control sites during 2018, suggesting that downstream effects of urban pollution were diminished after the WWTP upgrade. The assemblages at these impacted sites in 2018 were correlated inversely with concentrations of most algal and cyanobacterial biomarkers, as well as variation in river discharge (Fig. 5). Finally, bacterial assemblages in downstream sites on the Qu’Appelle River grouped together irrespective of years and were correlated positively with densities of most primary producers (cyanobacteria, green algae, and diatoms) and discharge. Taken together, planktonic bacterial communities were regulated mainly by NH₄⁺ pollution until the BNR upgrade, after which time more diverse regulation occurred due to variation in NO₃⁻, SRP, and in-stream conditions (phytoplankton, discharge, and turbidity).

Figure 5.

RDA of bacterial communities by site for each sampling year. Sampling locations are colour coded. Grouping of plot points towards the head of physico-chemical vectors (black arrows) indicates a positive correlation with the given parameter. Predictors are represented as arrows within RDA plots, whereas samples are plotted as points within ordination space. SRP is soluble reactive phosphorus. For 2012, the first axis explained 26.1% of the variance, while the second explained 20.1%. For 2018, the first axis explained 17.0% of the variance and the second explained 7.1% of the variance.

Bacterial community diversity and response to physical-chemical variables

Shannon Weiner estimates of bacterial diversity declined at wastewater-impacted sites after, but not before, the BNR upgrade (Fig. 6A), with GAMs explaining 52.1% of the deviance in temporal and spatial patterns. Similar trends were recorded for analyses of bacterial diversity based on species richness (49.9% deviance explained), Inverse Simpson (43.5%), and Chao1 indices (32.3%), although spatial patterns were less distinct with these metrics (Supplementary Fig. S22). Although mean diversity values were similar across years, patterns of bacterial diversity exhibited less spatial variation among adjacent sites in 2012 and instead showed a progressive increase from Wascana Creek headwaters to downstream locations in the Qu’Appelle River (Fig. 6A). In contrast, diversity patterns in 2018 revealed a marked increase upstream of the WWTP, a sharp decline downstream of effluent outfall, and another rise to consistent values within the Qu’Appelle River. These spatial patterns appeared to be consistent throughout May to September during 2018 but not in 2012, when diversity declined from a springtime maximum (Fig. 6D). Phytoplankton abundance declined during the 60 km interval downstream of WWTP (Sites 3–5) outfall in 2018, whereas abundance values increased progressively with distance downstream during 2012 (Fig. 6B). Phytoplankton and diversity had comparable patterns in each year, with an increase in both bacterial diversity and Chl a in 2012 from reference to downstream sites and more of a bimodal pattern for both bacterial diversity and Chl a in 2018, with peaks at reference and downstream sites.

Figure 6.

Spatial and temporal distribution of bacterial diversity (Shannon) and chlorophyll a in 2012 (blue) and 2018 (orange). Deviance explained for models was 52.1% (Shannon) and 77.8% (Chl a). Shannon diversity was selected for its greater deviance explained compared to other diversity indices (Richness, Chao1, and Inverse Simpson). Plotted lines are predicted mean values from ten thousand model simulations, and shaded regions are 95% credible intervals. Locations of sites are included on the x-axis in panels A and B. Note that site six is excluded from this model because this site is upstream of the confluence of Wascana Creek and the Qu’Appelle River and is therefore not part of the linear continuum.

Analysis of bacterial diversity in both 2012 and 2018 using GAMs confirmed that variation in bacterial diversity increased positively with discharge and abundances of diatoms and cyanobacteria, but declined with elevated pH, green algae abundance, and specific conductance (Fig. 7). When both years were included, the effects of exposure to NH₄⁺ were minimal until concentrations exceeded 1 mg N L⁻¹, after which diversity appeared to decline at higher N concentrations characteristic of tertiary-treated wastewater (>10 mg N-NH₄⁺ L⁻¹). Individually, NH₄⁺ was not a significant predictor in either year. Increased discharge favoured elevated bacterial diversity only during high flow regimes of 2012, whereas diversity decreased with specific conductance, pH, and green algae (Fig. 7B). During 2018, diversity was unrelated to low river flow, but increased with the abundance of diatoms, cyanobacteria, and secondarily chlorophytes (P = 0.1), while declining with specific conductance and pH (Fig. 7C). Diversity was controlled by a combination of physical-chemical variables (temperature, pH, turbidity, and discharge) and correlated to phytoplankton abundance.

Figure 7.

The relationship of water column bacteria diversity to physico-chemical and phytoplankton variables for all years (A), 2012 (B), and 2018 (C). Predictors include discharge, ammonium, pH, SRP, specific conductance, diatom abundance, chlorophyte abundance, and cyanobacterial abundance. Deviance explained for A, B, and C were 71.0%, 81.8%, and 76.1%, respectively. Plotted lines are predicted means and shaded regions are 95% credible intervals. The y-axis for the Shannon diversity are centred at zero. Flat lines with a slope of zero are not significant predictors.

Discussion

The objective of this study was to determine how the composition and diversity of planktonic bacteria varied in response to influxes of urban wastewater, responded to declines in NH₄⁺-N pollution following BNR upgrades, and were influenced by spatial and temporal gradients in phytoplankton abundance and composition, as well as in situ physical-chemical conditions. Before BNR, NH₄⁺ concentrations were elevated for 60 km downstream of the WWTP (sites 3–5), likely due to limited nutrient uptake and transformation because of nutrient saturation (Waiser et al. 2011), whereas NH₄⁺ was low after BNR (2018) and NO₃⁻ occurred at comparable levels in 2012. The overall composition of suspended bacterial communities at effluent-impacted sites became more similar to those seen in reference headwater reaches after the WWTP upgrade (Figs 2 and 5), suggesting recovery from the impacts of excessive NH₄⁺. Shannon Weiner diversity of planktonic bacteria declined at impacted sites relative to headwater reference locations in 2018 and in contrast to the linear increase from headwaters to downstream sites observed in 2012 (Fig. 6). The composition of bacterial genera shifted from those commonly found in high-nutrient environments (Betaproteobacteria, Rhodocyclaceae) to more oxic, nutrient-poor environments (Flavobacterium). Nitrifying bacteria (Nitrospira, Nitrosomonas) were difficult to detect at all sites, but were more regularly detected at all sites in 2012 compared to only impacted sites (3–5) in 2018 (Atashgahi et al. 2015). Interestingly, pathogenic taxa (e.g. Legionella) were observed downstream of wastewater inputs only after BNR. Despite this exception, our findings suggest that BNR reduced the deleterious effects of tertiary-treated urban effluent and helped shift the composition of bacterial consortia towards that seen in naturally productive or restored eutrophic streams (Riley and Dodds 2012, Bernhardt et al. 2018). However, these improvements were limited to a 60-km reach downstream of the WWTP (Sites 3–5), suggesting that other intensive activities (grain farming, livestock, hydrological management) may also influence bacterial composition in the downstream Qu’Appelle River (Hosseini et al. 2018, Swarbrick et al. 2022).

Wastewater influx and bacterial composition

Changes in bacterial community composition and, to a lesser degree, specific phyla and genera were consistent with river recovery at effluent-impacted sites following BNR. Specifically, ordination analysis revealed that bacterial community composition at impacted sites became more similar to those at upstream reference sites after BNR (Figs 2 and 6), suggesting a return to more natural stream microbiota (Philippot et al. 2021), even though the main phyla (Proteobacteria, Actinobacteria, and Bacteroidetes) were consistently abundant before and after the WWTP. Notably, these phyla are commonly enriched in urban effluent and gut flora (Wu et al. 2010, Atashgahi et al. 2015, Price et al. 2018, Martín and Liras 2020). There was also a more widespread occurrence of nitrifying bacteria in all sites in 2012 compared to 2018, suggesting that nitrification was not as important after BNR upgrade, which was observed in work by Dylla (2019). However, nitrifier levels in plankton bacteria in the water column were generally low and difficult to detect in both years; instead, nitrifier bacteria may had a greater role as part of periphytic or sediment assemblages, which is an avenue for future research (Drury et al. 2013, Sonthiphand et al. 2013) and would be expected in the Regina effluent. Specifically, the genera Flavobacterium (from the phylum Bacteroidetes), Mycobacterium (Actinobacteria), and Legionella (Proteobacteria) are all associated with pathogenic wastewater inputs (Ashbolt et al. 1995, Chahal et al. 2016), whereas Rhodocyclaceae (Proteobacteria) may have been associated with higher nutrient conditions before the upgrade (Atashgahi et al. 2015). Further, there was also a more widespread occurrence of nitrifying bacteria in all sites in 2012 compared to 2018, suggesting that nitrification was not as important after BNR upgrade despite the predominance of NO₃⁻, as seen elsewhere by Dylla (2019). However, nitrifiers were uncommon and difficult to detect in planktonic bacterial assemblages in both years. We speculate that nitrifying bacteria may have been more abundant in periphyton or sediment assemblages (Drury et al. 2013, Sonthiphand et al. 2013), although this hypothesis requires further study. Overall, changes in gross bacterial community structure, as well as those of specific taxa to some degree, were observed only in the effluent-impacted reach between the City of Regina above the confluence of Wascana Creek and the Qu’Appelle River, similar to patterns seen for suspended and benthic phototrophs in this fluvial system (Bergbusch et al. 2021a,b).

Ammonium concentrations appeared to be the primary control of bacterial community assemblage prior to BNR upgrade, whereas interactions with phytoplankton (diatoms, cyanobacteria) and physical-chemical conditions (pH, conductivity SRP, and NO₃⁻) were more influential drivers of community assembly after WWTP upgrade (Figs 5 and  7). Initially, we had expected that bacterial communities would be regulated principally by changes in algal and cyanobacterial densities, followed by the effects of physical-chemical conditions (Kent et al. 2007, Wang et al. 2015). In contrast, we found that nutrients were the paramount correlate of bacterial community composition until NH₄⁺ was removed (Peterson et al. 2011, Lee et al. 2017). After BNR, bacterial communities were correlated more with algal and cyanobacterial abundance, discharge, pH, and specific conductance, similar to controls seen in upstream reference locations and within the Qu’Appelle River (Figs 2, 6 and 7). In both field (Peterson et al. 2011) and lab settings (Grossart et al. 2007), changes in bacteria consortia have been shown to correlate with algae and cyanobacteria mainly under low-nutrient conditions, with exudates from phototrophs playing an important role in supplying nutrients to heterotrophic bacteria. However, bacterial communities after BNR in impacted reaches may also be influenced by bacterial communities suspended within the treated effluent (Atashgahi et al. 2015, Price et al. 2018), which make up almost all of the flow (30%–70% of the volume) of Wascana Creek after the spring freshet (Bergbusch et al. 2021a). This was not something we investigated directly and could benefit from further analysis.

Bacterial diversity and nitrogen removal

Bacterial diversity appears to have exhibited differential responses to wastewater influx depending on the year of study (Fig. 6). Elsewhere, responses of bacterial diversity to wastewater influx have also been inconsistent, with the influx of secondary treated effluent (low DOM, high NH₄⁺, and SRP) reducing (Drury et al. 2013) or increasing bacterial diversity (Wakelin et al. 2008). Here we observed both patterns, with relatively high diversity at impacted sites due to the influx of tertiary-treated effluent (low DOM and P, high NH₄⁺) and lower values after BNR was initiated (DOM, NH₄⁺, and P removed). Given that concentrations of NH₄⁺ were reduced >90% but that levels of NO₃⁻, dissolved P, and other parameters did not vary greatly due to WWTP upgrade (Table 1), we infer that bacterial diversity responded mainly to changes in the quantity of chemically reduced N. Elevated diversity during tertiary effluent treatment may reflect the sequential replacement of taxa capable of nitrifying NH₄⁺ by those that denitrify excess NO₃⁻ (Waiser et al. 2011, Dylla 2019), whereas the latter group only should be common after BNR. Future research with gene-specific primers will be required to validate this hypothesis and evaluate how the metabolic capabilities of bacterial assemblages change downstream (Peterson et al. 2011, Sonthiphand et al. 2013).

When considered across both years, absolute bacterial diversity was greatest in reaches with elevated river discharge and phytoplankton density (diatoms, cyanobacteria), but was correlated inversely with NH₄⁺ at high concentrations, pH, planktonic chlorophytes, and conductivity (Fig. 7A). Comparison between years showed that the strength of individual relationships varied through time and was likely related to changes in effluent solutes and their effects. For example, discharge was only a significant predictor of bacterial diversity during the high flow of 2012 (Fig. 7B), possibly reflecting increased influx of taxa or chemicals from land due to elevated runoff (Chen et al. 2018) or the scouring of benthic taxa and resuspension into the water column (Berger et al. 1996). Of interest, bacterial diversity increased in the Qu’Appelle River despite elevated cyanobacterial densities during summer 2018 (Bergbusch et al. 2021a,b), in contrast to Huang et al. (2022), who report that bacterial diversity declines due to cyanobacterial blooms in eutrophic Chinese rivers. The mechanisms underlying diversity correlations with pH and conductivity are less clear, however. Taken together, these findings suggest that N removal caused a decline in bacterial diversity between years, but that the many drivers of diversity (phytoplankton interactions, physical-chemical conditions) were important in both years.

Remaining risks of toxicity and pathogen exposure

Despite successful NH₄⁺ removal and changes to planktonic bacteria, BNR wastewater still poses a risk to downstream aquatic environments because of residual nitrogen and possible pathogen contamination. In our study, we found that BNR still allowed some NO₃⁻ release, albeit at a lower N:P ratio. Therefore, while not as severe, the effects of eutrophication may still persist in these P-rich waters, while a low N:P ratio could favour production of the cyanobacterial toxins (e.g. microcystin, anatoxin) by downstream cyanobacteria (Orihel et al. 2012, Hayes et al. 2019). We also observed potential pathogens in downstream Wascana Creek sites after the BNR upgrade, suggesting that some effluent-related microbiota remained after BNR or that there were additional, unmeasured contributions from livestock (Ashbolt et al. 1995, Chahal et al. 2016). Further investigation of the prevalence of potential pathogens in these streams is needed.

Conclusion

Fewer than 30% of cities employ BNR to remove N from urban effluent, even in developed countries. In part, slow uptake of this technology may reflect scientific and management concerns about the potential detrimental effects of selective N removal on eutrophication of downstream freshwaters, even though recent research demonstrates that NH₄⁺ removal reduces symptoms of eutrophication (excess primary production, toxic cyanobacteria) even in N-limited streams. Here we extend this observation and show that the removal of NH₄⁺ by BNR decreased bacterial OTU diversity but restored a composition similar to that seen in headwater reference sites. In particular, bacteria communities within 30–60 km of the WWTP outfall (Sites 3–5) exhibited a unique composition when polluted with tertiary-treated effluent but began to resemble headwater communities after BNR upgrades removed NH₄⁺ (Figs 2 and 5). Further changes in bacterial composition were correlated to a wider range of factors (NO₃⁻, SRP, and phytoplankton) rather than mainly NH₄⁺ concentrations, a pattern consistent with unpolluted eutrophic streams, or the return of lotic ecosystems to healthier conditions. While this study clearly demonstrates the effects of N pollution and its management on bacterial assemblages, future studies could use gene-specific probes to quantify changes in the metabolic capabilities of bacteria and better assess how the removal of NH₄⁺ via BNR impacts bacterial abundance, heterotrophic metabolism, and nutrient biogeochemistry.

Author contributions

Nathanel T. Bergbusch (Conceptualization, Data curation, Formal analysis, Methodology, Writing – original draft, Writing – review & editing), Alicia R. Wong (Formal analysis, Methodology, Writing – original draft), Jennifer N. Russell (Conceptualization, Formal analysis, Methodology, Writing – review & editing), Vanessa J. Swarbrick (Formal analysis, Methodology, Writing – review & editing), Claire Freeman (Conceptualization, Methodology, Writing – review & editing), Jordyn Bergsveinson (Formal analysis, Methodology, Writing – review & editing), Christopher K. Yost (Conceptualization, Funding acquisition, Methodology, Writing – review & editing), Simon C. Courtenay (Funding acquisition, Methodology, Supervision, Writing – review & editing), and Peter R. Leavitt (Conceptualization, Data curation, Formal analysis, Funding acquisition, Methodology, Resources, Supervision, Writing – review & editing)

Aknowledgements

We thank members of the Limnology Laboratory and D. Bateson for assistance with data collection and analyses. This study was conducted on Treaty 4 territory. We appreciate the willingness of the Indigenous Peoples of Saskatchewan to protect and share Canada’s water resources. This study is a contribution of the Qu’Appelle Valley long-term ecological research programme (QU-LTER).

Conflict of interest

The authors declare no conflicts of interest.

Funding

This work was supported by an NSERC post-graduate scholarship, an NSERC Alexander Graham Bell Canada graduate scholarship, NSERC Canada Discovery Grants programme, Canada Research Chairs, Canada Foundation for Innovation, the Province of Saskatchewan, and the University of Regina.

Data Availability

The data that support the findings of this study are available from the corresponding author on request.

Author notes

Present address: Environment and Climate Change Canada, Saskatoon, Saskatchewan S7N 3H5, Canada