Spatial and temporal patterns of urea content in a eutrophic stream continuum on the Northern Great Plains

Urea can degrade water quality and stimulate toxic phytoplankton in P-rich lakes, yet little is known of its sources, abundance, or transportation in lotic systems, particularly within the Northern Great Plains. We measured physico-chemical parameters biweekly during May–September 2010–2012 at 16 stations along a 250 km lotic continuum to quantify spatial and temporal variation in urea concentrations and discharge, and to identify potential regulatory processes. Urea concentrations were similar to those in regional prairie lakes (range 5.2–792.1, median 78.6 μg N L−1) with variable seasonal mean (± SD) concentrations (96.6  ± 96.1 μg N L−1) and fluxes (4.22 × 105 ± 257.6 μg N s−1). Landscape analysis with generalized additive models explained 68.3% of deviance in urea concentrations, with high temporal variability predicted mainly by positive relationships with nutrient content and chlorophyte abundance, but not temperature, dissolved organic matter, bacterial abundance, or urban effluent. Seasonal analysis revealed that during spring, urea content was correlated negatively with leguminous forage cover (% area) and positively with stream discharge, oilseed and cereal crops, and shrubs or deciduous plants, while during summer, urea concentrations were correlated negatively with discharge and leguminous crop cover, as well as nutrient levels. Mean porewater urea concentrations (528.5 ± 229.8 μg N L−1) were over five-fold greater than stream concentrations, suggesting that hyporheic production may offset declining influx from terrestrial sources during summer. We conclude that urea may be ubiquitous in eutrophic prairie streams and that management of its export from land may reduce detrimental effects on downstream lakes.


Introduction
Dissolved organic nitrogen (DON) is increasingly recognized for its role in lotic nitrogen (N) cycles and its impacts on water quality (Berman and Bronk 2003;Wiegner et al. 2006;Stanley and Maxted 2008) in lake (Berman 2001;Berman and Bronk 2003;Belisle et al. 2016) and coastal ecosystems (Seitzinger and Sanders 1997;Anderson et al. 2008;Korth et al. 2011). Of the diverse forms of DON, urea has received recent attention Solomon et al. 2010), as it is the most widely-used N-based fertilizer  and has been linked to the proliferation of harmful algal blooms (Bogard et al. 2020;Swarbrick et al. 2020). Indeed, experimental fertilization with urea increases phytoplankton abundance in phosphorus (P)-rich systems by 400% (Finlay et al. 2010;Donald et al. 2011;Bogard et al. 2020), and selectively stimulates toxic cyanobacteria (Solomon and Glibert 2008;Donald et al. 2010;Solomon et al. 2010). The effects of urea pollution are of particular concern in the Northern Great Plains (NGP), where the compound comprises [ 70% of N-based agricultural fertilizers (Statistics Canada 2016) and may be exported to P-rich downstream lakes (Bernot and Dodds 2005;Johnson et al. 2013;Rattan et al. 2016). Despite forecasts that urea supply rates from agriculture and urban sources will double by 2050 (Millennium Ecosystem Assessment 2005), little is known about the transport of urea from watersheds to downstream aquatic environments, particularly in the NGP.
While prior research on export of urea from watersheds to aquatic ecosystems has been limited mainly to temperate regions Han et al. 2015;Jackson 2016;Kibet et al. 2016;King et al. 2017), there are several reasons to expect that fluxes to freshwaters may be particularly elevated in the grassland environments of the NGP. First, these prairie watersheds include a variety of urea sources, including fertilizer applications (Davis et al. 2016;Kibet et al. 2016), intensive livestock operations (Peterson et al. 2004;Kibet et al. 2016), and urban effluents Bogard et al. 2012;Cozzi et al. 2014). Second, application rates of urea in the NGP are substantially elevated compared to other farming regions where urea export to streams has already been demonstrated (Belisle et al. 2016;Statistics Canada 2016). Third, the timing of urea application in NGP (spring, fall) is concurrent with hydrologic and meteorological conditions that may favor urea runoff (Stepanauskas et al. 2000;Siuda and Chróst 2006;Thorén 2007;Solomon et al. 2010). For example, spring snowmelt accounts for over 75% of annual hydrological runoff in Canadian Prairies , as well as up to * 80% of annual total nitrogen (TN) export to surface waters (Cade-Menum et al. 2013;Corriveau et al. 2013;Rattan et al. 2016). Further, cool temperatures in spring and fall can suppress enzymatic hydrolysis of urea to ammonium (NH 3 /NH 4 ? ) and favour export of undegraded molecules (Therkildsen and Lomstein 1994;King et al. 2017). Fourth, chemical inhibitors of urease are usually co-applied with urea to delay hydrolysis up to several weeks (Olson-Rutz et al. 2009), resulting in transport of undegraded urea to downstream ecosystems (Kibet et al. 2016). Together these findings suggest that streams of the NGP may exhibit particularly elevated concentrations and fluxes of urea during spring.
Urea concentrations within streams may also be regulated by in situ processes including production, transformation of nitrogenous compounds, and microbial uptake. For example, urea may arise from autochthonous sources including heterotrophic mineralization of particulate-and dissolved-organic N (PON, DON) within surficial sediments and water columns (Therkildsen and Lomstein 1994;L'Helguen et al. 2005;Thorén 2007; King et al. 2017), regeneration by suspended and benthic invertebrate grazers (Slawyk et al. 1998;Lomas et al. 2002;L'Helguen et al. 2005), photo-degradation of PON and DON (Bushaw et al. 1996;Bronk et al. 2010), and excretion from primary producers (Solomon et al. 2010). Such autochthonous urea production may be offset by heterotrophic and autotrophic microbial consumption, or by deposition as colloids or adsorbed onto particles. As well, the role of microbial consumption varies across environmental conditions, with potentially high rates of uptake by both heterotrophic microbes (Park et al. 1997;Jørgensen 2006) and autotrophs (Mitamura et al. 2000(Mitamura et al. , 2012Thorén 2007;Solomon et al. 2010;Donald et al. 2011;Bogard et al. 2020;Swarbrick et al. 2020). An integrated understanding of how these processes interact along lotic continua and across seasons is needed to assess the cumulative risks of urea pollution to surface water quality.
While dissolved inorganic nitrogen (DIN) export is well characterized (Mulholland 1992;Mulholland et al. 2000;Peterson et al. 2001;Seitzinger et al. 2002), less is known about the factors which regulate spatial and temporal patterns of DON and urea in lotic systems (Brookshire et al. 2005;Lutz et al. 2011;Johnson et al. 2013;Wymore et al. 2015;Pennino et al. 2016). Prior research shows that landscape patterns of DON are affected by regional and local land-use activities (Wiegner et al. 2006;Scott et al. 2007;Stanley and Maxted 2008), although few studies explicitly target transport of urea. Some conceptual models suggest that lotic ecosystems may act like an ''inert pipe'' when export is controlled primarily by hydrologic conditions (del Giorgio and Pace 2008). Under these conditions, urea concentrations may increase with catchment area and distance downstream. Alternately, lotic ecosystems may act as ''reactors'' (del Giorgio and Pace 2008; Kaushal et al. 2014), wherein metabolic processes are controlled by variation in biotic and physico-chemical conditions (Wiley et al. 1990;Young and Huryn 1996), as well as the degree of coupling of biogeochemical cycles of carbon (C), N and phosphorus (P) (Brookshire et al. 2005;Lutz et al. 2011;Wymore et al. 2015). In this case, concentrations of DON may either decrease if selectively consumed (e.g., 'Passive Carbon Vehicle Hypothesis'; Brookshire et al. 2005;Wymore et al. 2015) or may increase if DIN is preferred over DON by microbes (e.g., 'DON Release Hypothesis'; Lutz et al. 2011). Irrespective of the precise formulation, these models all suggest that urea export will be influenced by biological, physical and chemical gradients in both terrestrial and aquatic environments, as well as how these factors interact in time and space.
Here we collected stream samples every second week during May to September at 16 sites to quantify spatial and temporal patterns of urea concentration along a 250-km continuum of eutrophic prairie streams. Our study area spanned a gradient of agricultural and urban land use, enabling evaluation of the influence of both diffuse inputs from the watershed and point-source inputs from urban effluent. Our primary goal was to quantify the spatial (landscape) and temporal (seasonal) patterns of urea concentration and flux in grassland streams. Second, we sought to identify the relative importance of potential hydrologic, land use, and in-stream processes as controls of variation in urea concentration. We hypothesized that urea concentrations and fluxes would increase along the continuum (Heathwaite and Johnes 1996), but that concentrations would be low compared to those in regional lakes where sediments are a major source of urea to the water column (Bogard et al. 2012). Further, we hypothesized that urea concentrations would be highest during spring, when discharge and runoff is high (Martin and Harrison 2011;Corriveau et al. 2013;Rattan et al. 2016), temperatures are cool, and agricultural practices (fertilization, urease inhibitor application) increase urea concentrations on land (Swensen and Singh 1997;Siuda and Chróst 2006;Solomon et al. 2010). Finally, we hypothesized that receipt of urban effluent from a tertiary wastewater treatment plant (WWTP) would increase urea concentrations, particularly during late summer (Bronk et al. 2010;Bogard et al. 2012;Cozzi et al. 2014) when base flow is low and effluent is a large fraction of stream discharge (Bergbusch 2020).

Study area
The study area is centred within the Qu'Appelle River catchment, a region which comprises over 52,000 km 2 of mixed-grass prairie in southern Saskatchewan, Canada (Fig. 1). Wascana Creek is a first-order eutrophic stream that drains * 1400 km 2 of farmland before joining the main stem of the Qu'Appelle River. The headwaters of Wascana Creek are ephemeral, cease flow by mid-late summer, and empty into shallow Wascana Lake within the City of Regina. Below that point, Wascana Creek flows continuously, in part due to influx of tertiary-treated effluent (low P, elevated N) from the City of Regina (Waiser et al. 2011;Bergbusch 2020). Flow is also continuous in the Qu'Appelle River, which originates upstream of eutrophic Buffalo Pound Lake, and flows eastward to Pasqua Lake following its confluence with Wascana Creek (Fig. 1). The natural flow of the Qu'Appelle River is managed, with augmentation of discharge from Buffalo Pound Lake, and late-summer inputs from sub-saline Last Mountain Lake just east of the Qu'Appelle-Wascana confluence. Further details of system hydrology are provided by Haig et al. (2020).
Regional climate is characterized as cool-summer humid continental (Köppen Dfb classification), with short summers (ca. June-August; mean 19°C in July), cold winters (ca. November-February; mean -16°C in January), and low annual temperatures (* 2°C) with high seasonal variability (Leavitt et al. 2006;Haig et al. 2020). Average annual precipitation is * 380 mm of which * 30% falls as snow during the winter months ). Spring snowmelt, which typically occurs while soils are frozen during late March to early May, accounts for [ 70% of annual surface runoff (Fang et al. 2007;Dumanski et al. 2015). In contrast, rainfall occurs primarily in spring and early summer and decreases towards fall, causing low soil moisture and very low runoff by late summer ). In recent years, summer precipitation has accounted for an increased fraction of inflow to lakes in the Qu'Appelle River chain, in part due to declining snowfall (Vogt et al. 2018;Haig et al. 2020).
Land use within the study area is characterized by intensive agricultural development, wherein 95% of grassland has been converted to crop (dryland wheat, canola, mustard) or intensive livestock production (Hall et al. 1999). Wascana Creek receives stormwaters and wastewater effluent from the City of Regina (site 9 in Fig. 1; population * 210 500), as well as diffuse nutrient runoff from feedlot operations and golf courses. In many reaches, the riparian zones are degraded (Miki 2015), leading to designations of 'stressed' or 'impacted' for the lower Wascana Creek and Qu'Appelle River ecosystems (Davies and Hanley 2010).
Sixteen sampling sites were located along a * 250 km lotic continuum formed from Wascana Creek headwaters to the main stem of Qu'Appelle River near Pasqua Lake (Fig. 1). Seven sites were located upstream of Wascana Lake, five were located between Wascana Lake and the creek's confluence with the Qu'Appelle River, and four were located downstream on the Qu'Appelle River. One additional site (site 12) was located on the Qu'Appelle River upstream of the confluence with Wascana Creek. In addition to the 16 sampling sites, six gauged hydrometric stations operated by Water Survey of Canada (WSC), which recorded or modeled instantaneous discharge (m 3 s -1 ), were located within the study area ( Fig. 1). Field methods Ten study sites were sampled bi-weekly between May and September of 2010-2012 (Supporting Information), while an additional six sites were added in 2011 and 2012 ( Fig. 1; Supporting Information). Sites were sampled between 9:00 and 13:00 h. Samples were collected from each site as long as water was present, irrespective of discharge rates. Water temperature (°C), dissolved oxygen (mg O 2 L -1 ), conductivity (lS cm -1 ), and salinity (ppm) were measured in the middle of the water column at the stream thalweg using a YSI Model 85 m. Streamwater pH was measured at 5-cm depth using a calibrated handheld Oakton model 20 pH meter. A 10-L whole-water sample was collected from each site at mid-column depth along the thalweg. In 2011 and 2012, sampling was expanded to include turbidity and sediment porewater collection. Turbidity was calculated as the average of three replicate measurements using a portable LaMotte Model 2020we Turbidimeter. Standard volumes of sediment were collected for porewater analysis by inserting a polycarbonate tube into soft submerged substrate in * 0.25 m of water or at the deepest location. The proportion of substrate types was not recorded. After collection of a core with an undisturbed sediment-water interface, the overlying surface water was gently decanted, and the uppermost 15 cm of sediment transferred to a sterile Nasco Whirl-PakÒ, after which air was purged prior to being sealed. All water and sediment samples were kept on ice and returned to the laboratory for processing within two hours of collection.

Laboratory methods
Upon return to the laboratory, samples were filtered and preserved prior to analysis. Whole-water samples were filtered through GF/F glass fiber (nominal pore size 1.2 lm) and membrane filters (0.45-lm pore) before being frozen prior to nutrient analysis. Particulate organic matter (POM) on GF/F filters was frozen until analysis for chlorophyll a (Chl a) and biomarker pigments by high performance liquid chromatography (HPLC) using standard methods (Leavitt and Hodgson 2001). Subsamples were also refrigerated for analysis of dissolved inorganic-(DIC) and dissolved organic carbon content (DOC). Additionally, 10 mL of whole water was transferred using sterile technique to Vacutainer vials, preserved with 0.5 mL of Lugol's iodine solution, sealed and stored in the dark until analysis for microbial cell density the National Hydrology Research Institute (NHRI) in Saskatoon, Canada (see below). Porewater was extracted from sediment samples by centrifugation at 5000 rpm for 5 min (di Bonito et al. 2008), and supernatant was filtered through a 0.45-lm pore membrane filter, then frozen until analysis for nutrients.

Discharge rates and fluxes
Instantaneous discharge rate (Q) at the wadeable sites (sites 1-5; \ 1.3 m depth) was calculated by multiplying the measured water velocities within crosssectional increments and their area, following the twopoint and 6/10th depth methods of Buchanen and Somers (1976). Instantaneous Q at the non-wadeable sites was either recorded by provincial or federal WSA hydrometric real-time gauging stations ( Fig. 1), or calculated using the drainage area ratio method (Gianfagna et al. 2015). The drainage area ratio method involves adjusting the measured discharge rates from the nearest gauged station by an estimate of prorated discharge from the incremental effective drainage area (EDA) to the sampling site (Q incremental ), as well as any additional point source discharges or tributary inputs. The Q incremental is calculated by multiplying the discharge rate at the nearest gauged station by the ratio of the EDA of the ungauged sampling site to EDA of the gauged site as, where Q gauged is the measured instantaneous discharge at the gauged site, A incremental is the EDA between the sampling site and the gauged site, and A gauged is the EDA of the gauged site. Discharge rates calculated using the drainage area ratio method do not account for evaporative losses, groundwater inputs or losses, attenuation, or time of travel between gauged and ungauged sites, and are limited by the accuracy of measured EDA. Instantaneous Q was not calculated for site 14, as the site receives variable and ungauged hydrologic input from Last Mountain Lake ( Fig. 1). Although highly correlated (slope = 0.84, R 2 adj-= 0.856, P \ 0.0001), discharge measurements from wadable sites underestimated values calculated from gauging station data by about 0.75 m 3 s -1 at flows between 1 and 8 m 3 s -1 (Bergbusch 2020). Finally, instantaneous nutrient fluxes were calculated as the product of the instantaneous Q and nutrient concentration.
Watershed characteristics and geographic data processing The EDA of each sampling site was determined using the federal Agriculture and Agri-Food Canada (AAFC) Watersheds Project (2008) dataset. The AAFC dataset defines the EDA as the maximum area that could contribute runoff to an individual site during average hydrologic conditions (Mowchenko and Meid 1983), and provided delineated incremental EDAs for six of the study sites. Incremental EDAs for the remaining sites were delineated manually based upon the AAFC total EDA boundary layer, 1:50,000 topographic data, and 1:50,000 hydrologic data (Natural Resources Canada 2016) using ESRI ArcGIS 10.1. Because of the confounding effects of control structures that release or retain water from upstream lakes to the Qu'Appelle River, delineated EDA measurements for statistical analysis were restricted to sites 1 to 11. Annual land cover within the study area was extracted from the AAFC Annual Crop Inventory (2013), which uses optical (Landsat-5, AWiFS, DMC) and radar (RADARSAT-2) imagery to classify land cover into 65 different feature classes with a 56 m 2 (2010) or 30 m 2 (2011, 2012) spatial resolution. Annual land cover data within the total upstream EDA of each site were clipped from the AAFC dataset using ESRI ArcGIS 10.1. Once the individual feature classes within each EDA were extracted, they were aggregated into nine land cover classes (e.g., urban, pasture, crop type, natural vegetation, etc.).

Numerical analyses
Generalized additive models (GAMs; Wood et al. 2016) were used to estimate spatio-temporal trends in urea and its relationships to in-stream conditions. Models of urea concentration were fit with predictors including year, day of year (date hereafter), landscape position of the site (site hereafter), and in situ physicochemical parameters following Bergbusch (2020). Because GAMs capture both linear and non-linear relationships among variables (Pedersen et al. 2019), this analytical approach allows better estimation of trends in mean values and associated uncertainties. Herein range and null space of the smoothers' penalty matrices were fully penalized to perform variable selection, improve model fit and parsimony, and reduce potential concurvity of smooths (Marra and Wood 2011). We assessed basis size, dispersion of residuals, homogeneity of variance, and relationship between the observed and predicted response for all models to ensure model assumptions were not violated. Inference was based on tests for smooths and random effects developed by Wood (2013aWood ( , 2013b, plus the effective degrees of freedom of estimated smooth and model terms. Models were run in R (version 3.6.2; R Core Team 2019) using mgcv (version 3.6.2) with residual marginal likelihood (REML) smoothness selection (Wood 2011;Wood et al. 2017). We assumed urea concentrations were gamma distributed (positive, continuous responses).
Only variables that were collected in all three years were considered as predictors. Marginal smooths were visualized using ggplot 2 version 3.2.1 in R (Wickham 2016).
To additionally evaluate spatial and seasonal patterns of stream urea concentrations, data were aggregated by month and location. Sites were classified into one of three regions (Bergbusch 2020): headwaters, which comprised the reaches from the headwaters to the wastewater treatment facility (WWTP; sites 1-8); wastewater-impacted locations which included reaches from the WWTP to the confluence (sites 9-11), and; Qu'Appelle River sites, those immediately downstream of the confluence (sites 13-16). Urea concentrations were log 10 (x ? 1) transformed and tested for normality using a Shapiro-Wilk test. A twoway ANOVA was used to identify whether transformed urea concentrations differed between months or regions, and if there was an interaction between the effects of month and region on urea levels. The Durbin-Watson test was used to assess first-order autocorrelation and ensure that data met the assumptions of independence required for a two-way ANOVA, while the Levene's test was used to assess homoscedasticity. Paired t-tests on log 10 (x ? 1) concentrations of urea, NH 4 ? , NO 3 -, and TDN at sites 9 and 10 were conducted to assess if the discharge of urban effluent had a significant effect on instream N concentrations. Finally, because instantaneous discharge data were not normally distributed and data transformations to meet normality were unsuccessful, a non-parametric Kruskal-Wallis test with post hoc Dwass-Steel-Chritchlow-Fligner test for pairwise comparisons was used to identify significant differences in discharge rates between months.
We evaluated monthly relationships between discharge rate and urea concentration using Spearman rank correlation, where the sign of the coefficient (r s ) indicated whether hydrological flux appeared to enhance or dilute in-stream urea concentrations. The relationship between urea concentration and relative percent abundance of each land cover type was evaluated separately by month using Spearman rank correlation. The relative strength and direction of relationships between the selected instream water quality parameters and urea concentrations was then assessed using Spearman rank correlation. These tests were conducted in SYSTAT v. 13.0, except for Spearman rank correlations, which were conducted in TIBCO Spotfire v. 7.0.

Urea concentrations and discharge rates
Urea concentrations varied substantially within each sampling year at individual sites (Fig. 3), with a trend to increased levels during late summer in headwaters, and declines through the sampling period in several downstream Qu'Appelle River sites (Fig. 4a). Despite this pattern, mean urea concentrations were similar in all months (range 49.9-81.3 lg N L -1 ), with slightly elevated landscape values in mid-summer. Similarly, there were no significant differences (F (2, 290) = 1.597, p [ 0.05) in annual mean urea concentrations in headwater (66.1 lg N L -1 ), effluent-impacted (74.1 lg N L -1 ), and downstream Qu'Appelle River locations (61.2 lg N L -1 ). In contrast, two-way ANOVA revealed a significant interaction between site location and month (F (6, 290) = 8.183, p \ 0.001), consistent with the increase in urea concentration through summer in headwaters, but not at downstream sites (Fig. 4a).
Instantaneous discharge rates (Figs. 2b, 4b) were significantly different among sampling months (Kruskal-Wallis test statistic = 30.61, p \ 0.01), with pairwise comparisons and post hoc Dwass-Steel-Chrichlow-Fligner tests indicating that both May and June were significantly greater than other months, but  Fig. S2). Spearman rank correlations showed that although urea concentration and instantaneous discharge rate were not significantly related over the entire study period (Fig. 3; r s = -0.04, p = 0.46), they were significantly correlated within each month, though the direction of the relationship varied over the course of the sampling period (Supporting Information   Fig. S4). Specifically, urea concentrations increased with discharge rates during May (r s = 0.49, p \ 0.01) and June (r s = 0.36, p \ 0.01), but declined with increased flows during July (r s = -0.45, p \ 0.01) and August (r s = -0.33, p \ 0.01), even though discharge was greater in the latter two months. In situ correlates of urea concentration Analysis using GAMs showed that overall urea concentrations were correlated strongly (68.4% deviance explained) to variation in physico-chemical conditions, nutrient content, and autotrophic abundance (Fig. 5). Both year and sampling date were highly significant predictors, while sampling site was not an important variable influencing urea concentration. Urea concentration increased with year (not shown), and exhibited greatest levels during midsummer (Fig. 5a). Similarly, urea levels declined with specific conductance, showed complex relationships with discharge and DIC (Fig. 5a), but were uncorrelated to temperature, DOC, pH, and stream at the landscape scale. In general, urea levels increased significantly with concentrations of dissolved NH 4 ? and NO 3 -, but exhibited weak unimodal relationships (p * 0.1) with SRP concentrations and TDP:SRP ratios. Among biotic predictors (Fig. 5c), urea content increased strongly with chlorophytes (as Chl b) and secondarily cryptophytes (as alloxanthin), declined slightly with elevated densities of cyanobacteria (as echinenone) and total phytoplankton abundance (Chl a), and were unrelated to densities of heterotrophic bacteria and hyporheic concentrations of urea (not shown).

Spatial and temporal patterns in nutrient fluxes
Mean fluxes of urea increased from headwaters to downstream sites (Fig. 2b) and, despite large seasonal changes in urea concentrations (Fig. 4a), were Generalized additive model explained 68.3% of deviance in landscape patterns of urea, with significant effects of (top to bottom) year, sampling date, discharge, specific conductance, nitrate (NO 3 -), ammonium (NH 4 ? ), soluble reactive phosphorus (SRP), total dissolved N to SRP ratios (TDN: SRP), and abundance of diatoms (fucoxanthin), chlorophytes (Chl b), cyanobacteria (echinenone) and total suspended phototrophs (Chl a). Solid lines represent the mean effect, while shaded areas indicate 95% credible intervals. Non-significant predictor variables are not presented, but include temperature, dissolved organic carbon, and bacterial abundance. Note that the y-axis is centred at zero, whereas the X axis varies with predictor. Ticks on X-axis represent data observations relatively consistent across the sampling period (Fig. 4b). Mean (± SD) fluxes of urea (4.22 9 10 5 ± 257.6 lg N s -1 ), NH 4 ? (3.68 9 10 6 ± 257.6 lg N s -1 ), and NO 3 -(4.45 9 10 6 ± 257.6 lg N s -1 ) were lowest in the upper headwaters (sites 1-5), and increased downstream, following patterns of stream discharge. Mean urea flux remained relatively consistent along Wascana Creek headwaters (sites 1-5), but increased downstream of Wascana Lake and, in particular, below the wastewater outfall point (site 9), where export increased by 450% relative to upstream samples. Despite this pattern, the greatest urea fluxes were recorded immediately above (site 12) and below the Wascana Creek-Qu'Appelle River confluence (site 13; Fig. 2b). In comparison, the largest increase in mean NH 4 ? flux occurred immediately downstream of the wastewater treatment facility, with steadily declining values downstream of that point. Mean NO 3 fluxes also increased below urban effluent outfall, but remained elevated downstream, with the largest fluxes observed downstream of the Qu'Appelle River confluence (Fig. 2b). As shown elsewhere, these patterns are consistent with intense microbial nitrification downstream of the wastewater treatment plant (Dylla 2019).

Land use patterns
Pasture and crops consistently had the highest relative (%) cover within the EDA of each site (Fig. 6 a-c), whereas wheat (up to 66% cover), lentils (53%), and canola (41%) were the most common crops. The relative proportion of the aggregated land cover classes varied along the lotic continuum, with cereals, canola, and pasture dominating the headwater EDAs, whereas deciduous trees and shrubs increased in the downstream portions of the study area. In contrast, urban land was concentrated predominantly between sites 7 and 8, where it comprised 40% of incremental EDA. Relative proportion of land cover type also varied by year. For example, legumes dominated crop cover in 2010 (Fig. 6a), while a late spring and wet conditions in 2011 restricted crop production and increased the relative abundance of both surface water/wetlands (up to 32%) and pasture (up to 43%) within the watershed (Fig. 6b). Dry conditions in 2012 led to a greater abundance of cereal crops (Fig. 6c). In-stream urea concentrations were correlated significantly to five of the nine aggregate land cover types, though relationships varied by month (Table 1). Significant positive correlations between relative abundance of canola and urea concentrations were observed in May and June (r s = 0.31-0.52, p \ 0.05). June urea concentrations were correlated positively with cereals (r s = 0.32, p \ 0.01), exposed land (r s = 0.40, p \ 0.01), and shrub-deciduous cover (r s = 0.36, p \ 0.01), and negatively to pastures (r s = -0.30, p \ 0.05). Although cereal and shrub covers remained correlated positively to urea concentrations in July (r s = 0.42, p \ 0.01), legume cover also predicted urea concentrations in July and August (r s = -0.27 to -0.51, p \ 0.05).

Discussion
Urea comprises * 70% of nitrogenous fertilizer application in the Canadian Prairies Statistics Canada 2016), yet little is known of the patterns and controls of urea export to lotic ecosystems of the NGP. Unexpectedly, urea concentrations along the lotic continuum were elevated compared to those recorded in regional lakes (Finlay et al. 2009;Bogard et al. 2012). While urea fluxes remained relatively consistent through time at individual stations (Fig. 4b), two distinct longitudinal patterns in urea concentration emerged during the sampling period (Fig. 4a, Supporting Information Fig. S1), each mediated by interactions between hydrology, land use, and in situ stream processes. Consistent with our initial hypotheses, spring concentrations of urea increased along the continuum and were correlated positively to the relative abundance of N-intensive crops within the upstream catchment, although only when runoff ) and agricultural fertilization applications are typically elevated (Cade-Menun et al. 2013). In contrast, urea concentrations during summer decreased with distance downstream and were inversely related to hydrologic fluxes, possibly reflecting offsetting dilution of elevated hyporheic sources of urea (Supporting Information Figure S1), as seen elsewhere (Thorén 2007;King et al. 2017). However, contrary to our second hypothesis, effluent discharge had no effect on urea unlike other nitrogenous compounds (Waiser et al. 2011;Bergbusch 2020). Instead, urea concentrations during summer were correlated positively with dissolved N concentrations (NO 3 -, NH 4 ? ) and abundance of chlorophyte algae, taxa which benefit from pollution with urban wastewater (Bergbusch 2020;Bogard et al. 2020). Given these patterns, and considering the detrimental effects of urea on water quality, cyanobacterial abundance, and toxicity (Donald et al. 2011;Krausfeldt et al. 2019;Swarbrick et al. 2020), our results suggest that surface water management in the NGP needs to consider lotic inputs of urea as a potential factor degrading prairie streams and lakes.

Urea concentrations compared to other lotic systems
Despite the widespread use of urea fertilizers Statistics Canada 2016), there is little empirical data on spatial and temporal variation in urea within lotic systems in agricultural regions. Urea concentrations in this study (range 5.2-792.1 lg N Table 1 Spearman rank correlation coefficients (r s ) between instream urea concentration and relative abundance of land cover type with the total upstream effective drainage area.
Date range Canola/Mustards Cereals Exposed Grassland Legumes Pasture Shrub/Deciduous Urban Water Longitudinal patterns in urea along the lotic continuum Downstream patterns of urea concentration were the product of complex interactions between allochthonous and autochthonous inputs, biological processing, and in situ physico-chemical conditions. Although urea concentrations can increase downstream in large, coupled, lentic-lotic systems (Bogard et al. 2012), longitudinal patterns of urea recorded here varied among seasons, with urea concentrations increasing from first-order headwaters to downstream sites during spring and early summer, but decreasing along the continuum during mid-to late summer (Fig. 4), resulting in little spatial pattern at the annual scale (Fig. 5). As seen elsewhere for TN and DON (Martin and Harrison 2011;Corriveau et al. 2013;Rattan et al. 2017), lotic urea concentrations during spring appeared to be mediated mainly by high runoff that exported urea from terrestrial sources (Table 1). We hypothesize that agricultural fertilization in spring, often used in combination with urease inhibitors (Malhi et al. 2003), may combine with seasonally elevated runoff ) and lower temperatures to favour export of undegraded urea from agricultural fields or unimproved lots (Table 1) and that these inputs should increase as a function of the total drainage area. In contrast, the downstream decline in urea concentrations during summer (Fig. 4a) may reflect offsetting effects of elevated hyporheic influxes at all sites (Supporting Information Fig. S1), combined with progressive increases in stream flow (minimal in headwaters) that dilute in situ urea concentration (Fig. 4a) while maintaining overall high rates of export (Figs. 2b,4b).
Correlates of urea concentrations GAM analysis of landscape patterns of urea concentration suggest a paramount effect of seasonality, discharge, conductivity, in situ nutrient content, and chlorophyte biomass on water column levels of urea (Fig. 5). However, unlike most of the dissolved nutrient and biological predictors, variation in physico-chemical conditions appeared to have complex, non-linear relationships with in-stream urea content when analyzed over three years and 16 sites. As noted below, non-parametric analysis of individual months reveals pronounced differences in predictive models among seasons (Supporting Information Fig. S4) and suggests that some of the non-linearity in responses arises from changes in regulatory mechanisms among seasons.
Spring-Unlike previous studies that observed positive relationships between agricultural land use and DIN content in prairie surface waters (Corriveau et al. 2013), we found that urea concentrations were only correlated positively to agricultural land cover during spring and early summer (Table 1). Here, elevated urea levels occurred concomitant with seasonal maxima in runoff Pham et al. 2009;Fig. 3) and increased inputs of allochthonous nutrients, similar to patterns seen for vernal DON export elsewhere (Stepanauskas et al. 2000;Siuda and Chróst 2006;Martin and Harrison 2011;Tzilkowski et al. 2018), particularly in regions experiencing intensive agriculture (Pellerin et al. 2006;Johnson et al. 2009). In addition to vernal inputs, relationships between urea and in situ variables indicate that urea was subject to only minor processing during May and June. For example, the observed relationships between urea, DOC and NO 3 were most consistent with predictions of the ''DON release hypothesis'' (Lutz et al. 2011;Wymore et al. 2015), wherein DIN is preferred by microbes over DON when DOM pools are sufficient. Specifically, the DON release hypothesis predicts that DON concentrations will increase with distance downstream, such as seen here in spring (Fig. 4a).
Elevated urea concentrations in prairies streams during May and June coincide with high relative abundance of canola and other oilseed crops (Table 1), similar to temporal patterns of DON content seen elsewhere (Kemp and Dodds 2001;Pellerin et al. 2006;Johnson et al. 2009;Kaushal et al. 2014). Not only do these crops have elevated N demand that is met by fertilization (McKenzie et al. 2004), but, unlike cereal crops, oilseeds often have a ''memory effect'' on nutrient export (Heathwaite and Johnes 1996), wherein leaves and residues continue to release DON even after harvest. Further, the most common canola production practices involve broadcast application of urea, with urease inhibitors or polymer coatings, prior to or during peak hydrologic runoff in spring and fall (Olson-Rutz et al. 2009) or shortly before rainfall events or irrigation to prevent degradation to NH 3 and volatilization (Thiessen et al. 2005;Kibet et al. 2016). In contrast, inverse relationships between in-stream urea concentrations and the abundance of pasture land, which primarily comprises N 2 -fixing alfalfa and legume forages (Agriculture Canada 1982), are consistent with the low levels of N fertilizer applied to pastures. Given that global production of canola and oilseeds crops is projected to rise (Food and Agriculture Organization of the United Nations 2002), these results suggest that future management strategies should control urea export to protect downstream water quality.
Summer-Unlike prior studies of regional lake ecosystems (Bogard et al. 2012), summer urea concentrations did not increase markedly with drainage area. Instead, urea concentrations observed during July and August were inversely related to lotic discharge, but correlated positively to indicators of microbial abundance (e.g., chlorophytes), suggesting that urea levels were affected by a combination of in situ sources and biological transformations (Fig. 5). For example, longitudinal patterns showed that summer urea concentrations were highest in the headwaters, concurrent with elevated Chl a abundance, and declined downstream with increasing bacterial abundance (Supporting Information Fig. S4). These relationships are consistent with studies showing urea and DON uptake can be regulated by microbial activity (Cho et al. 1996;Park et al. 1997;Berman and Bronk 2003;Jørgensen 2006;Tzilkowski et al. 2018), and that DON can comprise 70-100% of TDN consumed by autotrophic and heterotrophic communities (Seitzinger and Sanders 1997;Wiegner et al. 2006). Additionally, low stream discharge during summer may enable longer contact between solutes and the streambed, thereby favouring both increased hyporheic loading and bacterial mineralization of urea (Therkildsen and Lomstein 1994;Park et al. 1997;Thorén 2007;King et al. 2017) relative to spring. The significant, but opposing, correlations of urea content with chlorophytes and cryptophytes (positive) and cyanobacteria (negative) (Fig. 5) also suggest that urea levels are regulated by the balance between production and consumption of urea by diverse photoautotrophs in lentic systems (Mitamura and Saijo 1980;Mitamura et al. 2000Mitamura et al. , 2012Bronk et al. 2007;Solomon et al. 2010;Bogard et al. 2020). Taken together, these patterns are consistent with a presumptive shift in regulatory mechanisms from mainly physical controls of terrestrial runoff in spring to metabolic controls related to stream production in summer. This latter case is consistent with the ''Passive C Vehicle'' model (Brookeshire et al. 2005;Wymore et al. 2015), wherein the low C:N conditions characteristic of agricultural catchments (Heinz et al. 2015) favour preferential uptake of DON over DIN (Ghosh and Leff 2013).
Some land-use practices may continue to influence lotic urea concentration during summer (Table 1). For example, legume crops were correlated negatively with urea concentrations during July and August (r s = -0.27 to -0.51, p \ 0.05), consistent with reduced fertilization of N 2 -fixing crops and their use in regional crop rotation to regenerate soil N content (Przednowek et al 2004). In contrast, cereal crops were positively related to urea concentrations in July, a period which coincides with the increased water requirements of the cereal flowering phase (Alberta Agriculture and Forestry 2011). Although rarer than in other agricultural regions, regional irrigation can be used to supplement late season water requirements for cereal crops in conjunction with N fertilization (Alberta Agricultural and Rural Development 2013), particularly within the Qu'Appelle River valley. Together, these results suggest that agricultural practices may influence urea concentrations summer-long, and that beneficial management practices and fertilizer management plans should consider the entire spring and summer period.

Hyporheic urea sources
This study demonstrated that urea concentrations were greatly elevated within the hyporheic zone (Supporting Information Fig. S1), consistent with measurements of organic nitrogen concentrations in the hyporheic zone of other lotic ecosystems (Triska et al. 1990;Lorite-Herrera et al. 2009;Tzilkowski et al. 2018), and mass-balance models which propose that sedimentary urea is a key influx to prairie lakes (Bogard et al. 2012). Although correlations between water-column and porewater urea concentrations were not significant at the annual scale (Fig. 5) and varied markedly among seasons, mean concentrations within the hyporheic zone were five-to tenfold greater than those in the water column at all sites (Supporting Information Fig. S1). Based on the strong positive relationship between urea concentrations and the abundance of chlorophyte algae that predominate in regional benthic habitats (Bergbusch 2020), as well as the large pool of sedimentary urea, we speculate that water-column urea may originate mainly from hyporheic sources during summer (Therkildsen and Lomstein 1994;King et al. 2017;Tzilkowski et al. 2018). Future research should investigate this linkage, as well as that between the presence of phytobenthos (e.g., chlorophytes), urea production, and microbial consumption.

Effect of urban wastewater inputs on urea concentration
Contrary to our initial hypothesis and in contrast to studies in other regions (Cozzi et al. 2014;Mitamura et al. 1994;Switzer et al. 2008), there was no marked effect of urban effluent influx on urea concentrations under most hydrologic conditions. However, our findings are consistent with those of Bogard et al. (2012) who found that urea concentrations within Regina effluent were usually \ 70 lg N L -1 , levels similar to those observed upstream of the wastewater treatment facility (97 ± 146 lg N L -1 ). We hypothesize that the low and similar concentrations of urea in urban effluent and headwater Wascana Creek may reflect the fact that the Regina WWTP used tertiary treatment processes in which DOM content was greatly reduced by microbial uptake during monthlong incubations in slow-flowing clarifier tanks. Instead, urban effluent may have exerted mainly indirect effects on urea content in eutrophic streams by altering the presence of nitrogenous chemicals and downstream biota (Bergbusch 2020) which subsequently affected urea concentrations (Fig. 5).

Implications for downstream water quality
This study suggests that land-use practices, in-stream conditions, and benthic processes interact to regulate the content and flux of urea through eutrophic streams of the NGP. However, the relative importance of these elements appeared to shift seasonally as hydrologic conditions changed, with the continuum acting as an ''inert pipe'', transporting increasing fluxes of urea downstream during the spring, but acting like a ''reactor'' during the summer (del Giorgio and Pace 2008). Overall, our findings suggest that there were few changes in mean monthly urea export to lakes during the study period, as elevated in situ production may have offset decreases in allochthonous sources and hydrological discharge. This finding contrasts extensive research which shows that most TN and DON is exported from prairie catchments during spring snowmelt (Corriveau et al. 2013;Martin and Harrison 2011;Rattan et al. 2016) and suggests that urea pollution may degrade downstream water quality in any season (Swarbrick et al. 2020).
Improved understanding of the processes regulating export of urea from agricultural catchments is critical for predicting its impacts on surface waters in the NGP. Previous mass-balance studies inferred that lotic ecosystems may not supply a high proportion of urea to lakes in the NGP (Bogard et al. 2012), whereas our study infers that flowing waters can provide an unexpectedly consistent supply of urea despite seasonal changes in hydrology, land-use practices, sources, and in situ processes. Given this observation, we suggest that development of nutrient objectives for urea may be an effective instrument to manage water quality (Chambers et al. 2011), particularly in the NGP where lotic N concentrations are elevated (Waiser et al. 2011), in situ NH 4 ? and NO 3 levels exceed federal chronic toxicity guidelines (CCME 2010(CCME , 2012, and P-rich lakes are degraded by N influx (Leavitt et al. 2006;Swarbrick et al. 2020). As agricultural production and fertilizer use is projected to double by 2050 (Millennium Ecosystem Assessment 2005), an increased focus on urea biogeochemistry may be essential to protect surface waters in the NGP. VJS's PhD thesis and does not necessarily represent the opinions or official position of Alberta Environment and Parks. This paper is a contribution of the Qu'Appelle Long-term Ecological Research program (QU-LTER). We gratefully acknowledge that this research took place on Treaty 4 territory, homelands of the Cree, Saulteaux, Dakota, Nakota, Lakota and Metis peoples.
Author contributions VJS and PRL designed the study: VJS lead data collection and performed most analyses; NTB completed GAM analysis;, VJS wrote the first version of the manuscript; All authors edited the manuscript and approved the final version.
Data availability The data that support the findings of this study are available from the corresponding author upon reasonable request.
Code availability R code used in generalized additive models available from NTB on request.

Declarations
Conflict of interest The authors declare no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.