Phytoplankton biomass is central to many concepts and theories in river ecology and is a core metric for evaluating eutrophication status, but the origin of the algae in suspension is commonly overlooked. River algae originate from 2 groups with contrasting life histories-benthic organisms that are occasionally washed into the plankton (tychoplankton) and true-plankton species capable of actively growing while drifting (Reynolds et al. 1994, Rojo et al. 1994, Reynolds 2006). Istvánovics and Honti (2011) also showed that self-sustaining phytoplankton in rivers can be meroplanktonic, i.e., species adapted to spend part of their life cycle in benthic environments (e.g., Stephanodiscus and Aula-coseira; Reynolds and Descy 1996). At any point along a river, the drifting phytoplankton community is a combination of recently entrained benthic algae responding to conditions over time near that location and planktonic algae growing in response to varying upstream conditions over time. Taxonomic characterization of the phytoplankton community can differentiate these spatial and temporal factors, but routine monitoring programs typically only measure chlorophyll a (Chl a) concentrations as a measure of algal production and flux.
Simple measures of benthic and suspended Chl a concentrations lack spatial or temporal context about algal origin, constraining the interpretation of changes in algal biomass to a benthic or planktonic domain. Either algal biomass at a place is assumed to result from the interaction of limiting factors (e.g., light, nutrients) over time (Dodds et al. 2002), or those factors are assumed to interact over space and time while algae travel downstream (Basu andPick 1995, Lohman andJones 1999). In rivers containing mixtures of true plankton, tychoplankton, and meroplankton strategies, interpretation of limiting factors on algal biomass export must account for variability over both space and time to recognize the dual origins of suspended algae. Assessments of trophic status based on Chl a concentrations (Bennett et al. 2017) might be refined with knowledge of how watershed and river characteristics influence the flux (biomass per time) and origin (benthic vs planktonic) of algae in transport. The current study is motivated by the need for extracting the maximal amount of information from commonly available measurements of river Chl a concentrations. This effort can lead to a better understanding of the origin and fluxes of river algae and to a more nuanced interpretation of the status and causes of river eutrophication. We previously tested an approach for estimating the origin of algae in suspension (benthic vs planktonic) and the size of these 2 algal pools by comparing high-frequency measurements of Chl a concentration with river discharge during storm events (Peipoch and Ensign 2022). The advantage of concentration-discharge (C-Q) relationships is that they only require measurements of Chl a and dis-charge at 1 location but enable inferences and integration of algal sources and fluxes from throughout the river ecosystem upstream (Speir et al. 2024).
Various metrics have been proposed and successfully employed to characterize and interpret C-Q relationships of solutes and particles over the last 2 decades (Butturini et al. 2006, Godsey et al. 2009, Lloyd et al. 2016). These metrics can also be used to assess Chl a dynamics. For instance, Godsey et al. (2009) proposed that the logarithmic slope of a C-Q relationship (b) broadly indicates whether lateral fluxes are limited by material supply (source limited; b < 0) or hydrological transport (transport limited; b > 0). Because the contribution of lateral inputs (e.g., subsurface flow paths, overland runoff) to Chl a fluxes is negligible, we hypothesize that source limitation (low; negative b values) indicate an insufficient supply of benthic (from resuspension) or planktonic (from highly concentrated pools) chlorophyll, causing dilution of Chl a concentrations at high flows (Fig. 1). Conversely, transport limitation (high; positive b values) will indicate a more continuous supply of Chl a throughout the high-flow period. Additional metrics that characterize nonlinear properties of C-Q responses are the flushing index (FI) and hysteresis index (HI). The FI is intrinsically related to b, but it focuses exclusively on the concentration response direction during the rising limb of the hydrograph (Vaughan et al. 2017). This focus on rising limb data makes the FI index particularly sensitive to the onset of elevated water velocity and river celerity and indicates either an initial dilution of phytoplankton concentrations (FI < 0) or an initial flush of benthic material (FI > 0; Fig. 1). The HI index reflects both the magnitude and rotation of C-Q responses, with clockwise rotation (HI > 0) indicating rapid source mobilization or the exhaustion of Chl a supply and counterclockwise rotation (HI < 0) representing slow or delayed mobilization over space or time (Fig. 1).
Here, we use C-Q relationships and high-frequency Chl a and discharge data from the National Ecological Observatory Network (NEON) to explore the origin and dynamics of fluxes of total suspended algae in rivers of different ecoregions. We also broadly characterize benthic and planktonic communities across NEON watersheds of different size, channel slope, and hydrologic regime to further explore how species composition and growth forms can explain Chl a C-Q patterns (Fig. 1). By doing so, we also revisit concepts and theory related to algal productivity in river ecosystems of varying size and condition.
Since 2017, the NEON program has provided an extensive network of automated instruments and field sampling at 34 different aquatic sites across the United States. NEON data provide a rare opportunity to combine high-frequency measurements from sensors with physical samples collected at multiple sites of different ecoregions. Such combination of data is frequently available only for 1 or few sites. NEON provides the necessary data to translate some of the methods and research previously applied to single sites to a large spatial and temporal scale. Our approach in this study was to obtain all the available and reliable data from NEON aquatic sites and to analyze high-frequency measurements of Chl a concentration and turbidity during high flows with a common analytical procedure. Then, we related these high-flow patterns to watershed conditions (e.g., land use, hydrologic regime) and algal communities (e.g., Chl a abundance, species diversity) to elucidate the controls on river phytoplankton communities and fluxes.
We performed all data compilation and analyses using R statistical software (version 4.1.2; R Project for Statistical Computing, Vienna, Austria). First, we downloaded all the available water-quality (NEON 2023e) and quantity (NEON 2023b) data for each NEON stream site (n 5 27) during the period from January 2017 to September 2022 with the loadByProduct function of the neonUtilities package (version 2.2.1; Lunch et al. 2023). Water-quality data included in situ sensor data of in vivo Chl a concentration and turbidity, available as 1-or 5-min instantaneous measurements. Water-quantity data included continuous measurements of stream discharge calculated from a stage-discharge rating curve and sensor-based measurements of water surface ele-vation. We downloaded only data complying with NEON's final quality flag (FinalQF 5 0). We discarded observations categorized as inaccurate (FinalQF 5 1) for any of the 3 parameters (chlorophyll, turbidity, and discharge). In addition, we also discarded NEON discharge data classified as anything but the highest quality category by Rhea et al. (2023) and substituted reconstructed composite discharge series data provided by Vlah et al. (2024) for those time periods. A total of 26 NEON stream sites had data that met these criteria. Then, we split water-quality and quantity datasets into separate files containing all the compiled data per site and year. For each annual file, we used the period.apply function of the xts package (version 0.13.1; Ryan and Ulrich 2023) to average 1-or 5-min data into 30-min intervals and then merged discharge and Chl a concentration data at a 30-min interval frequency into single files containing raw data for each NEON site and year that would be used for further C-Q analysis.
In addition to C-Q data, we also compiled all the available data for surface water chemistry (NEON 2023a), periphyton/phytoplankton properties (NEON 2023b), and algal community composition (NEON 2023d) for the same 2017 to 2022 period for each of the 26 NEON sites. Surface water chemistry included concentrations of dissolved SiO 2 (DSi; also known as silica), dissolved organic C (DOC), NH 4
1 , NO 3
, soluble reactive P (SRP), total dissolved N (TDN) and P (TDP), total N (TN) and P (TP), and total suspended solids (TSS). NEON staff collects water samples in the thalweg near the water-quality sensor set of wadeable streams and from a depth of 0.5 m in nonwadeable streams. At each site, 12 grab samples/y are collected at regular intervals during the sampling season and 14/y are collected on an irregular basis to capture major flow events. Samples are filtered within 3 h of collection, kept cold, and shipped overnight to external laboratories for analytical processing. A more detailed description of sample collection protocol and analytical procedures can be found in DP1.20093.001 (NEON 2023a). Algal properties include benthic and suspended Chl a abundance, total C, N, and P content, and d 13 C and d 15 N signatures. Algal community composition consists of the taxonomic identification and counts of microalgae from benthic (cells/m 2 ) and water-column (cells/ L) samples. In both cases, algal samples are collected 3Â/ y preferentially between spring and autumn at all NEON aquatic sites. At each sampling event, 3 to 5 benthic samples are collected among stream habitats (riffle/run/pool) for the 2 most dominant periphyton communities growing on natural substrates: cobbles (epilithon), sediments (epipelon), sand (epipsammon), woody debris (epixylon), or plant surfaces (epiphyton). One sample of seston or 3 samples of phytoplankton are collected from the water column of streams and rivers, respectively. Littoral benthic communities are also sampled in large rivers, providing 5 replicates of periphyton properties per site and date. A more detailed description of sample collection protocol and analytical procedures can be found in DP1.20163.001 (NEON 2023c). The analysis of taxonomy, nutrient content, and isotopic signatures of seston samples (water-column samples in wadeable streams) was discontinued in 2018, but Chl a abundance continued to be measured until the end of the study period.
In addition to NEON data, we obtained land-use cover data (Dewitz 2021) for most of the sites and aggregated land-use categories as the % cover of forested (deciduous, evergreen, mixed, and shrub), agricultural (cultivated crops and pasture), and urban (open to high development intensity) area in the watershed. We also obtained estimated stream slope and stream order for each watershed from DelVecchia et al. ( 2023), and we estimated the accumulated stream channel length at each site with several functions in the nhdplusTools package (version 1.3.0; Blodgett and Johnson 2023) in R. We examined linear relationships between land use, watershed area, and median values of streamwater chemistry parameters among the selected NEON sites (n 5 26) with simple linear correlation analyses. Land use and watershed area were log transformed to meet linear model assumptions. We used Shapiro-Wilk tests to assess the normality assumption and Pearson's correlation coefficient (r) to evaluate the correlation models' fit. We used linear regression analyses to evaluate the effects of storm intensity on C-Q metrics at each site. Storm intensity was log transformed prior to each regression analysis. We used Shapiro-Wilk tests to assess the normality assumption, Levene's test to assess the model residual variance, and the coefficient of determination (r 2 ) to evaluate the regression model fit.
For each site and y, we examined C-Q relationships by using a slightly modified workflow from Millar et al. (2022) designed to process subdaily C-Q data and analyze stormspecific hysteresis patterns. The workflow separates discharge into baseflow and stormflow components, smooths the stormflow component, and defines an amount by which discharge must exceed baseflow to trigger the start and end times of storm events over the entire time of the input file. Then, each storm event is identified with a unique number, the storm's peak discharge is identified, and C-Q data are normalized and separated into rising and falling limbs. We visually inspected each C-Q plot of the identified storm events and discarded those containing data gaps, abnormal Chl a values (e.g., 100-300 lg/L), or complex (e.g., multipeaked) hydrograph patterns. For the remaining events, calculations of HI and FI indices are computed as described in Lloyd et al. (2016) and Vaughan et al. (2017), respectively. Briefly, the HI at each discharge interval is determined by simple subtraction of corresponding interpolated values in the falling and rising limb so that HI represents the average of HI values, varies from 21 to 1, and its absolute magnitude is proportional to the concentration difference between the rising and falling limbs (i.e., area within the loop). The FI is effectively the standardized slope of Chl a concentration over time from the beginning to the peak of flow at each event. Total discharge (m 3 ) and duration of the storm (h) are also calculated for each storm event.
We defined storm intensity (m 3 /h) as the quotient between the accumulated water volume (m 3 ) flowing downstream throughout the entire duration (h) of each storm event (Millar et al. 2022). We scaled storm intensity values by watershed area (m 3 h 21 km 22 ) when comparing storm intensity among sties. In addition, we estimated 3 other metrics on an event-by-event basis: 1) the difference in Chl a concentration between baseflow and peak discharge (D[Chl a]); 2) the computation of b values-the linear slope between log-transformed discharge and concentration data (Godsey et al. 2009)-and associated p-values for each storm event; and 3) the ratio of the CVs of Chl a concentration and discharge (CV C ∶CV Q ), which can help in discerning between chemostasis or C-Q decoupling in C-Q relationships when b ≈ 0 (Thompson et al. 2011).
We examined differences in algal community composition among NEON sites at 2 different taxonomic levels. We first compared the relative proportion of algal divisions based on cell counts at each site and habitat (benthic vs planktonic). Then, we assessed community structure of benthic and suspended algae with nonmetric multidimensional scaling (NMDS; 500 random starts, 999 iterations) and Bray-Curtis dissimilarity matrices generated from genus relative abundance data. We conducted all community data analyses with taxonomic identification and cell count data collected 3Â/y at each site by NEON. For both benthic and suspended algae, we used the metaMDS function of the vegan package (version 2.6-4; Oksanen et al. 2022) to conduct the NMDS analysis, the envfit function to assess the p-values for correlations between each genus and NMDS axis, and the anosim function to assess differences in algal community similarity among different groups. In these analyses we did not include any benthic genera contributing <0.1% or planktonic genera contributing <0.01% to total algal community abundance (cell counts).
To further assess the contribution of benthic sources to algal fluxes during high-flow events, we used a heuristic approach based on a Chl a mass balance at the watershed scale. We first estimated total channel area for each NEON watershed within the continental United States, 3 of the 26 included in this study, by using channel length and drainage area for each river segment identified in the National Hydrography Dataset Plus (discover_nhdplus_id function of the nhdplusTools package). We estimated channel width for each river segment from drainage area (A) as w 5 A ab , with a and b being hydraulic geometry coefficients with mean and SD equal to 0.47 ± 0.13 (Singh 2003) and 0.8 ± 0.057, respectively (Galster 2007). Based on these coefficients' values and variation, and assuming they are normally distributed, we generated 1000 estimations of channel area for each watershed (n 5 23). We also generated a dataset of the same size (n 5 1000) for mean areal benthic Chl a (in mg/m 2 ), assuming a log-normal distribution, parametrized with the available NEON data at each site. Areal benthic Chl a represented the sum of all algal communities (epilithon, epipsammon, epixylon, epipelon, and epiphyton) growing on natural substrates in each site and on each date. We estimated data (n 5 1000) for total benthic Chl a at the watershed scale (in g of Chl a) by multiplying randomly selected values from both areal benthic Chl a and channel area estimations. Simulated values of total benthic Chl a at the watershed scale were randomly compared with 1000 simulated values of Chl a fluxes per storm event based on the available NEON data at each site and assuming a log-normal distribution, as with benthic Chl a. We note that this approach assumes that in vivo Chl a concentrations accurately reflect in vitro (extracted) concentrations, when, depending on the amount of nonphotochemical quenching, the ratio of in vivo Chl a fluorescence and extracted Chl a concentration can vary substantially (Cullen 1982). However, given that Chl a flux per storm event and total benthic Chl a per watershed span across 6 orders of magnitude, we believe this reductionist approach provides valuable comparisons between the benthic algal biomass prior to storm events and the amount of algal biomass washed out by these high-flow events across watersheds of different size and condition.
The 26 NEON sites included in this study are distributed across 16 ecoclimatic domains, hereafter referred to as ecoregions, spanning from the Taiga to the Neotropical region, with mean annual temperatures varying from 21 to 257C and an absolute range in mean annual precipitation of 2180 mm (Table 1). Most of the sites were located in 1 st -and 2 nd -order stream network segments, with nearly ⅓ of sites located in larger (up to 7 th -order) stream network segments. Watershed size spanned 1.1 to 47,110 km 2 , and channel slope ranged from 0.3 to 12.37 (Table 1). Forest vegetation cover varied greatly among the watersheds, ranging between 0.2 and 99.6% and reflecting the dominance of herbaceous communities and bare soils in some ecoregions. The watersheds of only 2 sites (Lewis Run and Arikaree River) comprised >50% agricultural lands, and the most urbanized site was Posey Creek in the Mid-Atlantic, with 22% of its watershed area occupied by urban development (Table 1). Overall, NEON streams included in this study largely represented low-order networks with minimal human impact but with varying geomorphological and hydrologic conditions, from montane channels with snowmelt-driven flow regimes (e.g., Blacktail Deer Creek) to low-gradient channels (e.g., Kings Creek) and intermittent streams (e.g., Arikaree River).
Over the study period, differences in streamwater chemistry were generally larger among ecoregions than within them. Exceptions were the Prairie Peninsula and Mid-Atlantic ecoregions, which include 2 streams with contrasting land-use influences (Tables 1, S1). Agricultural land cover was strongly correlated with median concentrations of TDN and TN (excluding Arikaree River; r 5 0.73 and 0.98, respectively; Table S3), but not with P, across the study sites. Median concentrations of DOC and TSS were not strongly related to land-use cover but were positively related to watershed area (r 5 0.68 and 0.66, respectively; Table S3). Differences across ecoregions showed that DSi concentrations, which are critical for diatom growth, were highest in the Southwest region and Central and Southern Plains, as well as in the 2 NEON streams from the Atlantic Neotropical region (Table S1). Concentrations of bioavailable P were slightly correlated to DSi (r 5 0.40; Table S3), were lowest in the Appalachian and Colorado Plateaus, and were consistently lower and less variable across sites than concentrations of bioavailable N (Table S1). Accordingly, 69% of the sites had potentially P-limiting conditions for algal growth, with molar DIN∶SRP ratios >16, whereas the remaining sites were well below this general threshold, indicating strong N-limiting conditions. In contrast, DSi ∶DIN ratios were consistently >1, indicating no Si-limiting conditions for diatom growth in any of the study sites.
Mean benthic algal biomass varied from 1.2 to 85 mg Chl a/m 2 among all sites and independently of watershed area (Fig. 2A). Suspended Chl a concentrations varied from 0.2 to 5.4 lg/L among NEON streams (Fig. 2B). At most sites, epilithic algae were more abundant than episammon and epixylon communities per unit area, except for Red Butte Creek and the Ozarks Complex sites, in which epipelon and epixylon abundance were the largest contributors to benthic algae, respectively. Epilithon and episammon N∶P ratios varied much less than inorganic N∶P ratios in stream water when compared across all sites (Fig. S1A, B), indicating stoichiometric homeostasis and reaffirming N-limiting conditions for algal growth in those sites with low DIN∶SRP ratios. Suspended algae, by contrast, showed more stoichiometric balance with the N∶P ratios in stream water (Fig. S1C).
Community composition of benthic algae differed among ecoregions (ANOSIM R 5 0.43, p < 0.001; Table S4, Figs 3A, S2A, B). We found strong evidence for an overall gradient of increasing abundance of epilithic algae from more septentrional (i.e., northern) ecoregions to the central and southern ecoregions, in which we found higher relative abundance of episammon and epixylon taxa, including those in littoral communities of large rivers (see NMDS2 in Fig. 3A). Some of the most dominant genera in benthic habitats-Cocconeis, Leptolyngbya, Nitzschia, Navicula, Phormidium, and Pseudanabaena-were also abundant and important for explaining variation in seston and plankton communities among ecoregions (Table S2, Figs 3A, B, S2C, D). Nonetheless, community composition of seston in small streams was different than phytoplankton communities of large rivers (Table S4). Benthic and suspended algal communities were both primarily composed of chrysophytes and cyanophytes, with chlorophytes only contributing to the planktonic communities of southern rivers (Figs S3, S4). The relative abundance of rhodophytes was higher (10-25% in cells/m 2 ) in N-limited stream sites of the Pacific Southwest and Southern Rockies (Fig. S3).
Chl a flux, mobilization, and sources Median Chl a fluxes per storm event varied across 4 orders of magnitude, from 0.01 g Chl a km 22 event 21 in the most intermittent stream to 52.9 g Chl a km 22 event 21 in the Black Warrior River of the Ozarks Complex biome (Fig. 4A). We also found large storm-specific Chl a fluxes in streams of the Southern Plains characterized by gentle channel slopes and high Chl a abundance during baseflow (Table 1, Figs 2A, 3A). These sites showed some of the largest increases in Chl a concentrations during storm events compared with other sites with storm events of similar intensity (Fig. 4B). In contrast, mountainous watersheds in the Pacific Northwest showed much lower D[Chl a] values but some of the most intense storm events per watershed area (Fig. 4B). At Martha Creek, despite limited benthic Chl a abundance (Fig. 2A), intense storms led to similar Chl a flux per storm as in the more productive low-gradient streams in the Southern Plains (Fig. 4A). Overall, western streams in the Rockies and Sierras with lower mean annual temperatures showed lower storm-specific Chl a fluxes compared with the rest of the NEON sites.
In general, streams within each biome had similar median HI values and variance (Fig. 5). Small, forested watersheds (<10 km 2 ) with limited benthic algae (Fig. 2A) showed some of the highest HI values throughout the study period, indicating rapid mobilization and exhaustion of Chl a , with 3 sampling events/y. There were 5 replicates for benthic algae and 1 replicate for suspended algae per event. Error bars represent SE of the mean. Numbers within parentheses inside a specific column show the total number of sampling events for that site (only shown for sites with <10 independent sampling events). Note log-scale axis for suspended Chl a concentration (B). NEON sites on the horizontal axis are ordered from smallest to largest watershed size. Complete names for each site can be found in Table 1. sources during storm events in these headwater streams (McRae, Martha, and LeConte creeks and Walker Branch; Fig. 5). For most NEON sites, the interquartile range of HI included positive and negative values, suggesting both rapid and lagged mobilization of Chl a sources, depending on the storm event. Of these sites, streams with snowmeltdominated hydrology or frequent spring flooding showed positive correlations between HI values and storm intensity, suggesting faster Chl a mobilization (HI > 0) during times of high algal growth and severe flooding disturbance (Table S3, S5). In contrast, other sites in the Neotropics, Mid-Atlantic, Prairie, and Ozarks showed no relationship between HI values and storm intensity, pointing to other factors that may be influencing variation in hysteresis loops across high-flow events (Fig. S5).
Most headwater streams with higher HI values (Fig. 5) had comparatively lower FI values than other watersheds of similar size (Fig. 6A). In some storm events, these small and forested streams, such as Walker Branch and Martha Creek, had negative FI values due to Chl a concentration peaking and falling sharply after discharge began rising. On occasion, McDiffett Creek in the Prairie Peninsula biome also showed this rapid flushing of algae well before the storm peak flow, generating negative FI and C-Q slopes (Fig. 6A, B). Beyond these occasional patterns resulting in negative FI values, FI was negatively correlated to watershed area and was consistently positive across sites (Fig. 6A), indicating a strong role of benthic resuspension processes. Similarly, C-Q slopes (b) were negatively correlated to watershed area, with predominantly negative b values in the Black Warrior and Tombigbee rivers of the Ozarks biome (Fig. 6B). These were the only 2 sites with consistently negative C-Q values, which, combined with their higher HI and lower FI values, suggest dilution effects of proximal sources of planktonic algae upstream. Indeed, the planktonic communities at the Black Warrior and Tombigbee sites differed from the others (Fig. 3B) and differed from the benthic communities in the same rivers (ANOSIM R 5 0.44 and 0.42, respectively, p ≤ 0.001; Table S4). Among all sites, nearly ½ of the storms (45%) had more variable Chl a concentrations than discharge (CV C ∶CV Q > 1), which suggests the existence of additional factors other than streamflow influencing streamwater Chl a concentrations. In particular, Chl a concentrations varied 2 to 3Â more than discharge in streams of Appalachian and Pacific Northwest ecoregions, in which benthic and suspended Chl a abundances were low (Fig. 2A, B) and HI indices were high (Fig. 5), leading to a positive correlation between HI and CV C ∶CV Q ratios across all sites (r 5 0.59; Table S3).
Finally, our heuristic approach to compare variation in storm-specific Chl a fluxes at each site with the total benthic Chl a supply in the watershed indicated an isometric scaling relationship (b 5 0.996; Fig. 7) between exported Chl a fluxes and benthic Chl a biomass that crossed 6 orders of magnitude for both fluxes and benthic standing stocks. This relationship between export and benthic supply was on or near the 1∶1 ratio line for many sites. These results indicate that estimated Chl a fluxes were regularly near or lower than the total amount of benthic algal biomass available for transport prior to storm events, suggesting a dominant benthic origin of Chl a fluxes among the 26 NEON sites in this study. Only the most arid watershed (Arikaree River) appeared to be a clear outlier, with much lower Chl a flux than that available in the watershed. This reduced flux was likely due to the intermittent flow regime in most parts of the watershed, which would result in a 1.
severe overestimation of benthic biomass by our heuristic approach.
Phytoplankton biomass at any point within a river network is shaped by proximal and distal factors affecting benthic and planktonic growth. This spatiotemporal complexity complicates relationships between river algal biomass and the factors governing its growth. Frequent flooding, rapid nutrient turnover, and the heterogenous light conditions of running waters magnify this complexity when measuring cause and effect, weakening stationary correlations between environmental conditions and algal biomass measurements (Chl a concentrations). As a step toward solving this measurement problem, our results showed that Chl a C-Q relationships efficiently characterized upstream variation in sources and supply of river algal fluxes across multiple ecoregions and scales of stream networks.
Benthic (tychoplankton and meroplankton) and planktonic contributions to river phytoplankton populations are spatially heterogenous and highly dynamic processes, which are difficult to measure. However, the nature of river ecosystems renders unique opportunities to reduce complex network-scale processes of this kind into more practical C-Q metrics that can capture the origin and mobility of solutes and particles (Klein 1984, Godsey et al. 2009). There are numerous examples of the use of C-Q relationships and metrics to describe dynamic processes of nutrient fluxes (Rose 2003, Kincaid et al. 2020), C fluxes (Butturini et al. 2006), and sediment fluxes (Lloyd et al. 2016). Of the suite of metrics that C-Q analysis typically provides, we found FI and b values particularly insightful to distinguish benthic and planktonic contributions to algal fluxes during storm events. Previous work relied solely on negative FI indices (i.e., dilution effects) to identify a dominance of planktonic sources (Istvánovics andHonti 2011, Peipoch andEnsign 2022), but our results showed that negative FI values can also occur in small (<10 km 2 ) forested watersheds, likely because of rapid mobilization by flood wave celerity of benthic algae directly upstream of the measurement site. To avoid misinterpretation, one must ensure that negative FI values are matched with negative b values that consider both rising and falling limb data. For example, the Black Warrior River was the only NEON site that had consistently negative b and FI values. This river is impounded by multiple locks and dams that create several reservoirs used for drinking water and hydroelectric power supply, as well as favorable conditions for planktonic algal growth (Honti et al. 2008, Dokulil 2014).
Unfortunately, the low number of storm events with available and reliable data (n 5 3) for the Flint and Lower Tombigbee rivers limits our assessment of true-plankton dominance in large river sites, but other evidence suggests they would likely display similar C-Q responses. For instance, recent work on lowland, nutrient-rich rivers comparable to Flint and Tombigbee river sites have shown clear dilution effects on Chl a concentrations by high flows (Bruns et al. 2022). Among this study's sites, planktonic taxa (e.g., Cyclotella) were more abundant in the 3 larger rivers than in seston communities of smaller streams. In fact, differences in phytoplankton composition were larger when comparing streams with rivers (i.e., by watershed size) than among ecoregions with distinct environmental conditions (see Fig. 3B). In stream sites, diatoms dominated both benthic and planktonic communities, which can indicate a large tychoplankton contribution to planktonic communities (Rojo et al. 1994). This pattern was weaker for rivers, which had lower similarity between benthic and planktonic communities than in streams and a larger presence of cyanobacteria in the water column.
Considering the differences between benthic and planktonic communities in NEON rivers and the prevalence of concentrating effects (FI > 0; b > 0) observed in low-and mid-order streams, true-plankton taxa appeared restricted to Chl a fluxes in channels of 5 th -order or above. However, more research is needed to strengthen this evidence, given the skewedness of NEON data to low-order streams. Shifts from benthic to planktonic dominance in sediment-water contact area (Gardner and Doyle 2018) or from benthic to water-column N removal processes (Wang et al. 2022) have also been identified in rivers of 5 th -order and above. This threshold at 5 th -order channels is consistent with the tendency toward a slope of 0 in chlorophyll C-Q relationships (Fig. 6B), but all >5 th -order rivers in our study are affected by the presence of impoundments. Others have suggested that the trend toward C-Q slopes of 0 from low-to highorder channels is a consequence of rivers acting as chemostats, where the large contribution of catchment processes and variability in small streams becomes overwhelmed by the influence of instream processes in wider and deeper channels (Creed et al. 2015). Given that the contribution of subsurface flow paths or overland runoff to Chl a fluxes is negligible, we interpret values of b ≈ 0 (or slightly negative b values) as a consequence of similar dilution effects of baseflow plankton biomass and concentration effects by benthic resuspension and tributaries on Chl a fluxes. Inherently, this interpretation implies that although planktonic contributions increase with stream order, the contribution of benthic sources to phytoplankton biomass may still be substantial in >5 th -order free-flowing rivers.
Hysteretic responses of Chl a fluxes measured at a single location emerge from network-scale variation in source, availability, and mobilization processes of algal biomass pools upstream from the point of measurement. Sources are broadly characterized between benthic and planktonic; availability is proportional to biomass growth and abundance; and mobilization depends on near-bed velocities, shear stress, and lateral connectivity outside the channel. All these elements interact within the channel to shape C-Q responses of Chl a fluxes and are, in turn, affected by biogeochemical conditions that depend on climatic and land-use settings at regional to continental scales. For instance, light limitation is highly heterogeneous in small, forested streams (e.g., because of canopy openings and channel orientation; Warren et al. 2016, Burrows et al. 2021), minimal in midorder rivers (Julian et al. 2008), and again, highly variable in the water column of large and deep channels (e.g., because of sediment loads; Gardner et al. 2020). At the same time, light availability in river channels is ultimately constrained by latitudinal variation in sunlight exposure and regional gradients of forest cover. It is therefore necessary to consider both fine-and large-scale controls on river Chl a fluxes and their hysteretic patterns to elucidate the role of algal sources, availability, and mobilization.
Volume 44 June 2025 | 153
Hysteresis responses are intrinsically variable (Butturini et al. 2008), but we found reasonable similarity in Chl a fluxes and hysteretic responses between streams of the same biome (see Figs 3A, 4A). This within-biome similarity suggests that regional controls on river algal growth and dynamics are at least partially reflected by the timing and intensity of Chl a fluxes within each biome. Likely broadscale controls are resource availability (nutrients and light), hydrologic regimes (flooding intensity and seasonality), and sediment transport regimes (frequency and extent of coarse bed mobilization). In this study, positive effects of favorable nutrient and light conditions on algal biomass pools manifested as generally higher benthic and suspended Chl a abundance and D[Chl a] values in temperate ecoregions with low-density forest cover (e.g., Prairie Peninsula or Southern Plains) and in nutrient-rich ecoregions (e.g., Mid-Atlantic). On the other hand, ecoregions with more-intense spring flooding events led to signs of exhaustion of Chl a supply (positive HI) during this season. Chl a fluxes per storm were maximized when both abundant algal biomass and exhaustive mobilization coincided, such as in the Pringle Creek or Blue River sites.
However, storm intensity and Chl a abundance could also individually explain the quantity of Chl a flushed out of a watershed. We did not anticipate that Pacific Northwestern streams could have large Chl a fluxes per watershed area and storm compared with other ecoregions with higher light and nutrient availability. Benthic algal growth in these mountainous streams with dense canopy cover is strongly colimited by low light and nutrient availability (Warren et al. 2017), yet storm events in Martha Creek had comparable Chl a fluxes to streams in the Prairies and Southern Plains with higher benthic Chl a abundance. We speculate that storm events in small Northwestern streams mobilized most of the benthic algal biomass accumulated in the channel because of intense flooding and a tectonically active geomorphology that generates sufficient power to frequently and extensively mobilize and transport bedload, resulting in higher entrainment of benthic algal biomass. Their nearly 1∶1 match in our watershed-scale Chl a mass balance (Fig. 7) supports this hypothesis. In contrast, more productive lowland streams require forceful spring floods to mobilize all available algae, as suggested by the appearance of positive HI values only when storms were most intense.
Appalachian streams flushed out comparatively less Chl a per storm than Pacific Northwest streams despite their similar size, HI indices, and algal biomass abundance. A potential explanation is the higher resistance to near-bed shear stress by nonfilamentous diatom communities compared with filamentous algal communities (Biggs and Thomsen 1995). Nonfilamentous diatoms of Fragilaria and Cocconeis genera were more abundant in Appalachian streams, whereas filamentous cyanobacteria (Homoeothrix and Phormidium) were more dominant in the Northwest (Figs 3A,S3,S4). High presence of Cocconeis, Synedra, and Fragilaria diatoms has been associated with early successional stages of periphyton communities (Davie et al. 2012), which should be common in Appalachian streams where forest phenology provides very limited time and opportunities for algal growth. Therefore, while clockwise Chl a hysteresis suggests rapid mobilization of benthic algal pools in both ecoregions (Fig. 1), stronger periphyton resistance to near-bed velocities in Appalachian streams could explain the reduced Chl a fluxes. In other cases, however, variation in HI indices and Chl a fluxes between ecoregions was not consistent with differences in community composition. For instance, benthic algal communities were similar in streams with N-limiting conditions, including the Pacific Northwest streams, but Chl a fluxes and HI indices in Pacific Northwestern headwaters were substantially larger than in other N-limited watersheds. It is possible that the high intensity storm events in the tectonically active Pacific Northwest are responsible for the different HI and flux patterns in these streams.
Chl a fluxes are primarily controlled by storm intensity, but local effects and spatial heterogeneity in periphyton abundance, community composition, or scouring resistance can also affect algal biomass fluxes. We generally found higher variability in Chl a concentrations relative to discharge (i.e., CV C ∶CV Q > 1) in small watersheds with limited algal abundance, more positive C-Q slopes, and clockwise hysteresis and less variability in larger watersheds. Overall, variance in Chl a concentrations decreased faster than discharge variance with watershed area, in contrast with higher CV C ∶CV Q ratios in larger watersheds typically observed for dissolved constituents (Diamond andCohen 2018, Marinos et al. 2020). For solutes, high CV C ∶CV Q ratios combined with positive C-Q slopes can suggest the existence of dynamic processes connecting sources and streamflow over short time scales, whereas high CV C ∶CV Q ratios and near-zero C-Q slopes may indicate instream processing or transformations during the high-flow event (Speir et al. 2024). For algae, we interpret that high CV C ∶CV Q ratios and C-Q slopes indicate intense and variable inputs of benthic algal biomass being scoured by high flows, whereas high CV C ∶CV Q ratios combined with near-zero or even negative C-Q slopes may arise from benthic resuspension and planktonic dilution processes interacting over short time scales.
We found an isometric relationship between benthic standing stocks and fluxes of Chl a, with Chl a fluxes being equal to roughly 18% of the estimated benthic Chl a in the river network (a 5 0.18; Fig. 7). This proportion is reasonable, considering some of the assumptions in our simple mass-balance model. First, we compared Chl a fluxes with an entire channel network area serving as a potential source of Chl a supply, regardless of the size of the contributing watershed or the storm event. Second, a complete removal of all benthic algae in each site would be necessary for a site to appear on the 1∶1 line of the relationship, but high flows can remove between 10 and 80% of algal biomass, depending on water velocity (Flinders and Hart 2009). For instance, Bott and Newbold (2023) documented average removal of benthic Chl a between 69 and 81%, depending on storm intensity, in a 3 rd -order stream. Accordingly, our initial expectation was that Chl a fluxes would be much lower than estimations of benthic abundance, given that single storm events may neither affect all of the watershed nor remove all benthic chlorophyll. However, we expected these discrepancies in the mass balance to increase substantially with watershed size because of a larger overestimation of benthic Chl a abundance and potential removal. Instead, the observed isometric scaling suggests a constant ratio of benthic to exported Chl a across a large gradient in channel area.
A potential explanation for the consistent discrepancies in mass balance across watershed size is that the isometric relationship simply reflects isometric scaling between discharge and watershed area, factors that directly affect Chl a fluxes and benthic abundance, respectively. This explanation would add to other evidence in our study already suggesting a dominant role of benthic contributions to phytoplankton biomass in that the amount of algal biomass flushed out is consistently proportional to in-channel benthic abundance, regardless of the size of the contributing watershed. Larger watersheds experience larger floods, but if the relative influence of flood disturbance on the removal of benthic biomass remains constant across channel size, that would leave algal biomass fluxes as proportional to the amount of benthic material that can be mobilized. Recent work modeling the autotrophic biomass recovery post-flooding disturbances found no relation between the relative magnitude of flow disturbance and watershed size or order (Blaszczak et al. 2023), showing faster recovery in wider channels with limited flow regulation favoring shallower waters, higher light availability, and faster algal growth (Lowman et al. 2024). The assumption of a relatively constant proportion of benthic biomass removal across different stream sizes is consistent with the Blaszczak et al. (2023) finding of faster algal biomass recovery in river channels that offer better conditions for algal growth and storm resilience.
Alternatively, isometric scaling between benthic and exported Chl a could emerge from an increasing contribution of planktonic and off-channel algal sources in larger watersheds, which were not included in our mass-balance model. These unaccounted sources are likely more relevant in large and impounded river-floodplain systems, where empirical work has shown allometric scaling of dischargearea relationships (Galster 2007), suggesting other mechanisms behind the observed isometric relationship than a purely hydrologic scaling phenomenon. It is plausible that phytoplankton biomass accumulation prior to storms and algal production in off-channel habitats can compensate for the larger gap we expected to find between in-channel benthic abundance and Chl a fluxes in high-order watersheds. The contribution of in-channel phytoplankton biomass to measured Chl a fluxes is negligible in small-order watersheds but potentially relevant in larger ones, such as the Blue, Flint, and Black Warrior river sites. A simple way to assess the potential contribution of baseflow phytoplankton biomass to Chl a fluxes during storms is to assume a mean baseflow water depth of 1 m for all siteswhich likely over-and under-estimates the actual depth in lower-and higher-order watersheds, respectively-and to compare benthic and suspended Chl a concentrations in Fig. 2B. For smaller streams, areal benthic Chl a would be equivalent to 10 to 30Â higher suspended Chl a concentrations, indicating that the unaccounted in-channel suspended biomass can only make a negligible contribution to Chl a fluxes flushed out during storms. For larger watersheds, even with a conservative 1-m depth estimation, benthic and suspended Chl a concentrations would be quite similar. Moreover, several studies have supported the importance of natural in-channel storage zones (e.g., backwaters) hosting much higher phytoplankton concentrations than the river's main flow area, with strong implications for phytoplankton recruitment (Reynolds 1995, Reynolds andDescy 1996). Much less is known about the contribution of off-channel habitats in the floodplain (e.g., side channels, springs, wetlands, or ponds) to longitudinal changes in river phytoplankton concentrations, but the mechanisms ought to be similar to those of in-channel storage zones when high flows connect them to the river's main flow.
In conclusion, our results show that Chl a C-Q responses provide valuable information on the origin, abundance, and mobility of river algae. Spatial and temporal changes in the location and abundance of algae in a watershed are a fundamental discrepancy between the river continuum (Vannote et al. 1980) and patch dynamics concepts (Townsend 1989). We contend that a series of highfrequency Chl a sensors deployed in nested locations within a large watershed could help elucidate the interaction between local and routing controls on algal dynamics proposed by these theoretical concepts. In general, we conclude that consistent patterns in Chl a C-Q responses suggest a ubiquitous dominance of benthic supply to phytoplankton biomass. Pending further analysis of C-Q data from large, free-flowing watersheds or nested designs, we also conclude that phytoplankton dynamics in rivers largely resemble that of streams, unless flow is severely regulated (as in the Black Warrior River). Phytoplankton fluxes in a particular reach seem to mainly persist by continuous input of benthic supply (tychoplanktonic and meroplanktonic strategies), a process that is obviously accentuated during high-flow periods but also likely dominant during baseflow conditions. Future quantitative assessments of the benthic dependence of river phytoplankton would contribute to better understanding the role of algal production and algaegrazer interactions on sustaining higher trophic levels of >5 th -order rivers.
Author contributions: MP: conceptualization, methods, data analysis, writing-original draft; MD and SHE: conceptualization, methods, writing-review and editing.
This research was financially supported by the National Science Foundation grant DEB-2213574 to MP, MD, and SHE. Comments and suggestions from 2 anonymous reviewers and the journal's technical editor substantially improved this manuscript.
| Origin and fluxes of river algae M. Peipoch et al.
| Origin and fluxes of river algae M. Peipoch et al.
Volume 44June 2025 | 149
| Origin and fluxes of river algae M. Peipoch et al.
Volume 44June 2025 | 151
| Origin and fluxes of river algae M. Peipoch et al.
| Origin and fluxes of river algae M. Peipoch et al.
| Origin and fluxes of river algae M. Peipoch et al.