Nutrients influence seasonal metabolic patterns and total productivity of Arctic streams

The seasonality of gross primary production (GPP) in streams is driven by multiple physical and chemical factors, yet incident light is often thought to be most important. In Arctic tundra streams, however, light is available in saturating amounts throughout the summer, but sharp declines in nutrient supply during the terrestrial growing season may constrain aquatic productivity. Given the opposing seasonality of these drivers, we hypothesized that “shoulder seasons”—spring and autumn—represent critical time windows when light and nutrients align to optimize rates of stream productivity in the Arctic. To test this, we measured annual patterns of GPP and biofilm accumulation in eight streams in Arctic Sweden. We found that the aquatic growing season length differed by 4 months across streams and was determined largely by the timing of ice‐off in spring. During the growing season, temporal variability in GPP for nitrogen (N) poor streams was correlated with inorganic N concentration, while in more N‐rich streams GPP was instead linked to changes in phosphorus and light. Annual GPP varied ninefold among streams and was enhanced by N availability, the length of ice‐free period, and low flood frequency. Finally, network scale estimates of GPP highlight the overall significance of the shoulder seasons, which accounted for 48% of annual productivity. We suggest that the timing of ice off and nutrient supply from land interact to regulate the annual metabolic regimes of nutrient poor, Arctic streams, leading to unexpected peaks in productivity that are offset from the terrestrial growing season.

Gross primary production (GPP) in streams and rivers represents a critical source of energy to food webs (Brett et al. 2017) and serves as an important driver of biogeochemical cycles (Heffernan and Cohen 2010;Rocher-Ros et al. 2020). The significance and fate of GPP are partly determined by its seasonal timing and predictability, which can guide patterns of community composition and consumer growth (Tonkin et al. 2017;Huryn and Benstead 2019), while also creating time windows when nutrient uptake and retention are optimized . At network scales, distinct patterns of GPP among streams also create emergent productivity patterns (Koenig et al. 2019) that have consequences for mobile aquatic consumers (e.g., Saunders et al. 2018), as well as for ecological interactions across individual reaches (e.g., Uno and Power 2015). While we know much about what controls instantaneous rates of GPP (e.g., Hill et al. 2009), we are only starting to explore how multiple physical, chemical, and biotic factors interact to shape seasonal and annual patterns (Bernhardt et al. 2018).
Traditionally, ecologists have viewed incident light (e.g., Finlay 2011) and hydrologic disturbance (e.g., Uehlinger 2006) as primary factors driving temporal change in GPP. At network scales, light regimes are regulated by channel width and orientation (Vannote et al. 1980;Julian et al. 2008), but can be further modified by the seasonal progression of riparian canopy cover (Hill and Dimick 2002). Differences in light regimes can drive distinct seasonality in GPP, where larger, unshaded systems assume patterns that mirror the terrestrial growing season , while smaller, forested streams show predictable peaks in spring, prior to leaf out and canopy closure Lupon et al. 2016). Hydrologic disturbance, in the form of scouring flows or drought, may further modify the relationships between light and photosynthesis (e.g., Beaulieu et al. 2013), creating either stochastic (Blaszczak et al. 2019) or predictable (Ulseth et al. 2018) seasonal patterns in GPP. Together, light and disturbance regimes can serve as effective organizers of metabolic regimes among stream and rivers ). Yet, efforts to resolve the drivers of GPP seasonality have so far been applied to a geographically limited set of systems (e.g., Appling et al. 2018b)-as we apply such analyses more broadly, the emergence of additional drivers seems likely.
A third candidate driver of seasonal productivity is the temporal variation in nutrient availability to primary producers. Indeed, seasonal variation in terrestrial nutrient supply to streams is widespread (Stoddard 1994;Sponseller et al. 2014), as is evidence for nutrient limitation of benthic algal growth (Keck and Lepori 2012) and GPP (Mulholland et al. 2001;Bott et al. 2006). Despite these observations, nutrient supply has not been incorporated into seasonal GPP models, in part because technological and/or logistical constraints often prevent pairing nutrient data with continuous metabolism records (Ulseth et al. 2018;Mejia et al. 2019). Additionally, many long-term metabolic records come from streams/rivers that are to some degree influenced by land use activities (e.g., Appling et al. 2018a) and may be sufficiently enriched with nutrients that they have little or no effect on productivity patterns. However, for the vast, remote regions of Earth drained by nutrient poor waters, temporal variability in terrestrial supply of nitrogen (N) and phosphorus (P) may act as an overlooked driver of GPP seasonality.
Nutrient supply might be a particularly important regulator of GPP seasonality in Arctic tundra streams. Given ice/snow cover throughout the dark Arctic winter, light limitation sets the coarse boundary condition for stream growing season lengths (Huryn et al. 2014). However, throughout spring and summer, sparse canopy cover likely result in light saturation of benthic surfaces (Hill et al. 2009). Superimposed on this light regime is a strongly seasonal pattern of stream nutrient concentrations, with relatively high levels in winter, and lowest during summer (e.g., Holmes et al. 2000;Mann et al. 2012), potentially due to catchment nutrient retention during the terrestrial growing season (Weintraub and Schimel 2005b;Ohte 2012). As a result, tundra streams are often strongly nutrient limited when incident light is greatest, as evidenced from studies in Alaska (Peterson et al. 1983), Northern Sweden (Myrstener et al. 2018), Iceland (Friberg et al. 2009), andGreenland (Pastor et al. 2020). At the same time, elevated nutrient concentrations in spring and autumn may create time windows that support elevated GPP. Overall, while a few studies report estimates of GPP from Arctic streams (e.g., Huryn et al. 2014;Song et al. 2018), the lack of continuous annual records limits our understanding of productivity regimes at such high latitudes. This knowledge gap is all the more problematic given the ongoing climatic changes in the Arctic related to longer growing seasons (ACIA 2004), modified landwater exchange of resources (Frey et al. 2007;Kendrick et al. 2018), and altered hydrological regimes (McClelland et al. 2006).
Here we explore the seasonal and annual patterns of GPP in an Arctic river network. We first ask what factors regulate the timing and magnitude of GPP throughout the year. In addition, we assess how seasonal patterns of GPP relate to differences in total, annual productivity. Finally, we evaluated how multiple productivity patterns among streams lead to an emergent pattern at the network scale. To answer these questions, we measured daily GPP and monthly periphyton biomass together with light, temperature, nutrients, and discharge at eight streams within the Miellajokka drainage network in Northernmost Sweden. We hypothesized that the shoulder seasons-the transitional periods between winter and summer-are critical time windows when light, disturbance, and nutrient availability align to optimize rates of stream productivity in these systems.

Study site
Miellajokka is a 52 km 2 catchment located 200 km north of the Arctic Circle in Sweden (Fig. 1). This area is dominated by Fennoscandian highland tundra, which is floristically similar to tundra ecosystems in other parts of the Arctic (Virtanen et al. 2016). While various terms are applied to this region (e.g., Arctic, sub-Arctic, Oroarctic), it lies within the Arctic boundaries defined by the Arctic Monitoring and Assessment Program on the basis of latitude, elevation, vegetation, and occurrence of permafrost (AMAP 1998). Annual temperature in this area averaged −1 C in the last climatic period  but is now above 0 C (Kohler et al. 2006). Annual precipitation is 350 mm yr −1 at the closest weather station, with about 60% falling outside the summer months (Kohler et al. 2006). However, the winter of 2017/2018 only had 70 mm of snow (climate data from Swedish Meteorological and Hydrological Institute [SMHI]). In this region, discontinuous permafrost exists at elevations above 800 m a.s.l, but the precise extent of this cover in the Miellajokka is unclear (Gisnås et al. 2017). The vegetation below 700 m.a.s.l, is dominated by mountain birch forest (Betula pubescens spp. Czerepanovii) with herbaceous understory and patches of shrubs and wetlands. Above the tree line, the catchment is characterized by tundra vegetation, including dwarf shrubs such as Empetrum hermaphroditum and Betula nana. Cryoturbated soils are typically found between 800 to 1000 m.a.s.l., and the top of the catchment (above 1200 m a.s.l) is bare moraine with permanent snow fields.
The hydrological patterns in this catchment are described by Lyon et al. (2018) and are shaped by a clear mountain-to-valley transition (e.g., Covino and McGlynn 2007). In short, above the tree line, the Miellajokka catchment is steeper and represents a water collection zone where the streams are dominated by precipitation inputs. At lower elevation, the catchment flattens and glacial deposits create an alluvial valley which supports diffuse flow through deep sediments. There, separate headwaters emerge, which are not connected to surface water from the tundra streams. These headwaters have comparatively stable hydrological and thermal properties, and finer benthic substrate composition. Discharging groundwater from the alluvial valley ultimately combines with stream water coming from tundra in the larger birch forest streams at the catchment outlet. We focused this work on eight sites that capture this variation in catchment structure and hydrology. Four are located in the tundra (M10, M9, M6, and M17), of which M17 is considerably smaller in drainage area (see Table 1). M18 and M16 are two small birch forest streams that emerge within the alluvial valley bottom. Finally M1 (the main stem) and M2 represents those larger birch forest streams near the base of the catchment that integrate hydrologic inputs from the tundra and alluvial valley.

Physical and chemical parameters
We recorded light and temperature at each site (HOBO model UA-002-64, Onset Computer Corporation) and converted lux to photosynthetically active radiation (PAR, μmol m −2 s −1 ) by a factor 0.0185 (Thimijan and Heins 1983). Light is presented as daily photon flux (DPF, mol m −2 d −1 ) using an average of two sensors at each site. We measured dissolved nutrients, including nitrate (NO 3 ), ammonium (NH 4 ), soluble molybdate reactive phosphate (SRP) every 5 weeks at all sites. Dissolved inorganic nitrogen (DIN) was calculated by summing NO 3 and NH 4 and is the N form used in all statistics and figures. On average, NO 3 accounted for ca. 80% of the DIN pool. Samples for nutrients were filtered in the field (0.45-μm Millex HA filter, Merck Millipore) and frozen before colorimetric analysis of NO 3 -N (ISO 13395:1996; Method G-384-08 Rev. 2), NH 4 -N (ISO 11732:2005; Method G-171-96 Rev. 12), and SRP (ISO 6878:2004;Method G-297-03 Rev. 1) with a QuAAtro39 (Seal Analytical). For the main stem (M1) we also included additional nutrient samples from spring and summer of 2017 collected through the Swedish Infrastructure for Ecosystem Science (SITES) program. These data were used to augment our assessment of seasonal nutrient patterns as well as the GPP-nutrient relationships at this site.
At each site, we had continuous measurements of water level (HOBO model U20-001-04, Onset Computer Corporation) at 30-min intervals. Site-specific water level was related to reach-specific depth from direct measures of depth at 5-10 transects. To obtain daily discharge (Q), 4-21 manual discharge measurements were conducted between 2015 and 2018 at various flow levels using the salt slug method (Moore 2005). Rating curves between reach depth and discharge were established for each site using linear models (R 2 between 0.83 and 0.99). As a measure of hydrological disturbance, we used the number of flood events above the 75th percentile of flow (Blaszczak et al. 2019) but excluded events in the smallest stream, M17 because of very low maximum flows (< 20 L s −1 ). Stream slopes were calculated based on a 2-m digital elevation model.

Biomass accumulation and standing stocks
We estimated biomass accrual of periphyton from chlorophyll a (Chl a) on ceramic tiles between August 2017 and September 2018. Five unglazed ceramic tiles (each 94 cm 2 ) were deployed every 5 weeks and retrieved after 5 weeks in each site throughout the year (except under dry winter conditions). At each site and sampling point we also estimated biofilm standing stocks from five cobbles. We excluded bryophytes from this because they were not clearly present at all sites. However, there was one stream (M16) with substantial bryophyte cover, which is captured by metabolic estimates but not by the biomass assessment. Cobbles and tiles were transported to the laboratory before biofilm was removed with a hard-plastic brush. The suspension of each individual cobble and tile was filtered through a nontreated GF/C filter (Whatman) for Chl a analysis. We estimated the area scrubbed on each cobble using aluminum foil and an area to weight relationship. Filters for pigment analysis were frozen at − 80 C before extraction. Chl a concentrations were analyzed spectrophotometrically following methods described by Steinman et al. (2017) including acetone extraction and acidification correction for pheophytins.

Metabolism modeling
Metabolism was estimated using the single-station diel oxygen method approach from August 2017 to September 2018 in all stream reaches. At M1 we were able to include additional metabolic estimates from spring and summer of 2017 using dissolved oxygen (DO) data from the SITES. In addition, we compiled daily estimates for M18 and M6 during summer 2016, which were only added to explore interannual variability but were not included in any statistical analysis.
At each site, we measured DO at 10-min intervals with optical oxygen loggers (miniDOT; Precision Measurement Engineering Inc.). We chose reaches that did not gain or lose more than 10% water over 50-500 m stretches to minimize groundwater and tributary inputs. However, differences in lateral inflows from day to day could still bias estimates of metabolism, although this is less of a problem for gross primary production than for respiration (Hall and Tank 2005). Each oxygen sensor was intercalibrated three times: before, during, and after deployment with DO saturated water (100% DO   Table 1 after air bubbling) and DO depleted water (0% DO using dry yeast) and there was no drift above instrumental accuracy (± 5%). GPP was estimated using Bayesian inverse modeling (Hall and Hotchkiss 2017) based on time series of DO, water temperature, light (from Abisko Scientific Research station), as well as daily estimates of gas exchange rate coefficient (K) and stream depth (z). The main equation used was The change in oxygen over time (g O 2 m −3 ) equals all oxygen produced by photosynthesis (GPP, g O 2 m −2 d −1 ) minus all oxygen consumed by respiration of both autotrophs and heterotrophs (ER, g O 2 m −2 d −1 ), and the rate of gas exchange between the water and air (K, d −1 ). We modeled three parameters (GPP, ER, and K) but enforced large constraints on K to minimize problems of equifinality (Appling et al. 2018a). To do so, we informed the model with K estimates and uncertainty boundaries that were based on prior propane releases and nighttime regression analysis (following Rocher-Ros et al. 2020). Briefly, we selected days when the nighttime regression produced K values with an R 2 > 0.6, and then constructed linear models between K and discharge to obtain daily estimates of K. At one site (M2) the relationship between K and discharge was weak (R 2 = 0.09) and so we used K values directly obtained from nighttime regression and linearly interpolated missing days. Additionally, the K vs. discharge relationship was rarely linear at the lowest discharge ranges, which typically occur under ice cover or close to ice cover. For these time periods, priors for K were manually set after visual inspection of model fits, aiming for parameter estimates that produced nonnegative GPPs for the last and first visible DO curves of the autumn and spring, respectively. Furthermore, the standard deviation (SD) of the prior distribution of K was increased from 0.05 SD to 0.1 SD for those time periods to allow a better model fit. Finally, we filtered data for erroneous model days by using the mean average error between the observed and the modeled DO concentrations. All days with a mean average error larger than 0.2 were removed (Lupon et al. 2019). When GPP is very low, a poor model fit can still produce a small error, thus all remaining days were also visually inspected to further exclude evidently erroneous estimates. Following these guidelines, we removed on average 21% of analyzed days across all sites, of which a large fraction were days during ice cover (see Fig. 2).

Data analysis
To aid in data analysis and interpretation, we distinguished specific time windows within which to explore patterns of GPP. First, we defined winter as the time when stream bottoms receive light < 1 mol m −2 d −1 , roughly corresponding to November to April. We then defined the open water season (May to October) as the period when most streams are ice free and receive light. Summer was defined as June to August, based on the terrestrial growing season in the area (Christensen et al. 2012). Finally, the shoulder seasons were defined as the transitional periods between summer and winter (May, September, and October, Fig. S1).
We used simple linear regression to compare the relative importance of temperature and locally measured light on daily GPP during transitions from autumn to winter and winter to spring. To organize this, we isolated the time windows when light decreased (autumn) or increased (spring) between 0 and 7 mol m −2 d −1 , which roughly corresponds to October and November in autumn and February to March in spring. We used the first and last day of detectable GPP as the starting and end points, respectively. For the winter-spring transition, only two streams gradually changed from 0 to 7 mol m −2 d −1 (M2 and M16) and could thus be analyzed in this way. Other streams in the network transitioned from no light to levels well above 7 mol m −2 d −1 in a matter of a few days at ice off and were not included in the analysis. All regression models were checked for possible serial autocorrelation in the residuals using the Durbin-Watson test in R (package "lmtest"). When this test was significant, we applied an autoregressive transformation (AR1) that adjusts linear models based on the strength of the residual autocorrelation (Cochrane and Orcutt 1949). We then used an additional Durbin-Watson test to ensure that the lag one residuals in transformed models were no longer autocorrelated. These transformed models have adjusted R 2 and p values accordingly.
We also used linear regression to assess the controls over changes in GPP throughout the open water season at each site. To do this, we calculated average rates of GPP for the 10 days around each grab sample for water chemistry and used this value in regression analyses against concentrations of DIN and SRP. Importantly, there was no autocorrelation between nutrients and other physical or chemical variables, except for a weak, inverse relationship between light and DIN at M1 (R 2 = 0.4, p = 0.04). Furthermore, Durban-Watson statistics indicated no temporal autocorrelation in the GPP-nutrient regression given the 5 week sampling lag in this test (p > 0.05), and we therefore applied no adjustments to these regressions. To explore how these bivariate relationships were organized across the network, we plotted the output from regression analyses (R 2 ) against the average DIN concentration during the open water season for each site. We chose DIN concentration as the organizer because previous studies in this region have documented N limitation for similar streams (Myrstener et al. 2018). Finally, for comparison, we used regression analysis to assess relationships between GPP vs. light and temperature throughout the open water season from daily data. In this case, we accounted for temporal autocorrelation in the regressions as described above.
We evaluated the drivers of annual GPP and biomass accrual (Chl a on tiles) during the open water season using partial least square regressions (PLSR) with the R package "pls" and leave-one-out cross-validation. PLSR identifies the relationship between predictor (X) and response variables (Y) through a linear, multivariate model and produces latent variables that maximize the explained variability in Y and reduce the original multidimensionality (similar to a principal component analysis). We opted for PLSR due to the high number of predictor variables compared to the low number of observations (Carrascal et al. 2009). The most important predictor variables were identified based on VIP scores (variables important in projection) above 0.9, according to Mehmood et al. (2012). For the GPP model, we included SRP, DIN, temperature, the number of flood events, winter season length, stream slope, and drainage area. In the biomass model, we included, SRP, DIN, temperature, light, and days post flood.
An additional goal of our analysis was to make an initial estimate of seasonal metabolic patterns at the network scale, accounting for differences in the size and areal coverage of our focal study streams. To do this upscaling, we made several simplifying assumptions. First, the total stream area was based on channel widths, measured manually at 168 locations along all major drainages in the network, and the total channel length was derived from a 2-m digital elevation model (see Rocher-Ros et al. 2019). To scale daily GPP from each oxygen logger to the upstream reach, we assumed that rates would not change dramatically upstream, unless catchment land cover changed from birch forest to tundra. Accordingly, we extrapolated GPP numbers linearly to the area of each reach for the smaller subcatchments (M18, M16, M17, M9, M6, and M10) that drained a single land cover type (i.e., tundra or birch forest). For the two larger streams (M1 and M2), which drain a mix of land covers, we applied measured rates of GPP for the upstream channel until land cover transitioned from birch forest to tundra. Above this transition point, and for any other tundra streams not captured by empirical estimates (only two, see Fig. 1), we assigned a daily GPP value calculated as the average across the three larger tundra streams (M6, M9, and M10), which we believe best represent this landscape type.
Finally, to compare these network scale patterns with the seasonality of terrestrial GPP, we compiled data from the NASA Soil Moisture Active Passive mission which produces GPP data at a 9-km resolution globally. GPP estimates are based on satellite measurements of surface micrometeorology, the fraction of absorbed PAR, surface moisture, and land cover class among other variables. The global data were downloaded using the package "smapr" for R (version 0.2.1), and we used an average daily GPP for 2017-2018 based on the three pixels that the Miellajokka catchment covers.

Stream physical and chemical variables
All streams displayed clear and similar seasonal patterns in light with levels below 1 mol m −2 d −1 by mid-November and the highest light levels during May and June (20-150 mol m −2 d −1 , see the seasonal light pattern of M1 in Fig. S2). However, the onset of light varied among streams due to differences in snow and ice cover (black bars, Fig. 2). For example, M16 reached 1 mol m −2 d −1 already in early February, whereas M17 did not reach this light level until the end of May. Overall, temperature averaged 7.5 C across streams during the open water season. During summer, the small birch forest streams (M16, M18) tended to be coldest (max = 9.5 C), followed by the larger mixed cover streams and smaller tundra streams (M1, M2, M6, and M17, max = 15.7 C), and then finally the larger tundra streams reached the highest temperatures (M9 and M10 max = 19.5 C).
Average summer discharge varied among streams from 10 (M17) to 1408 (M1) L s −1 , with the average of the other six being 90 L s −1 (SD ± 43). We observed that all tundra streams (M6, M9, M10, and M17), and one of the birch forest streams (M2) were dry at some point between January and late April. Additionally, there were two distinct spring flood events, the first in mid-May and the second at the end of June (see the seasonal discharge pattern of M1 in Fig. S2), but only the second flood peak propagated through all streams. Apart from the spring flood, the streams had a total of zero (M16 + M17) to four (M1 + M10) additional flood events during summer, where discharge increased above the 75th percentile. Most of the flood events corresponded to > 50% increase in discharge compared to the previous 6 days.
All streams showed similar seasonal patterns in DIN with fivefold higher concentrations in winter compared to summer considering all sites (average 147 vs. 31 μg N L −1 , Fig. S3A). At the onset of the terrestrial growing season there was a distinct drop in DIN at all sites except for the two small birch forest streams (M16 and M18) where concentrations remained elevated for another month. Average summer DIN varied among streams, with the lowest being 9 μg N L −1 (in M17 and M2) and the highest 53-68 μg N L −1 (in M16 and M18, respectively), while the other 4 streams averaged 29 (± 6) μg N L −1 . We had insufficient observations for a robust assessment of discharge-concentration relationships at all sites, but there was evidence for dilution of DIN at the main stem (M1), and to a lesser degree at M18, but not for the remaining sites (Fig. S4). Even though M1 showed signs of dilution, based on the larger set monitoring data at M1, the decline in DIN during snowmelt and throughout the summer was greater than would be expected only from dilution (e.g., when compared to more conservative solutes; Fig. S5). One consequence of this is that the low DIN concentrations during summer also corresponded to low fluxes at most sites. Finally, SRP showed little seasonal variation and was frequently below 2 μg L −1 at all sites (Fig. S3B), yet most sites did show a small peak in SRP during the spring flood.

Biomass accrual and standing stocks
Algal biomass accrual on tiles during 5-week deployments ranged from 0.01 to 5.3 μg Chl a cm −2 among all sites (Fig. S6A). Overall, the highest rates of accrual were observed in autumn, from September to October (average 0.92 ± 1.3 μg Chl a cm −2 ) and the lowest during June and July (average 0.13 ± 0.1. μg Chl a cm −2 ). Across sites, average accrual during the summer (June to August) was 0.34 μg Chl a cm −2 (± 0.8). Finally, Chl a standing stock on cobbles varied across sites and over time from 0.02 to 5.0 μg Chl a cm −2 , with similar seasonal patterns as observed for biomass accrual (Fig. S6B).
Overall, birch forest streams (M1, M2, M16, and M18) had higher biomass accrual compared to tundra streams (open water average; 0.8 vs. 0.2 μg Chl a cm −2 ). PLSR results suggest that this variation during the open water season was best explained by DIN concentration (correlation score = 0.81) on the first component, which explained 30% of total variance (Fig. S7). The second axis, driven by water temperature (correlation score = 0.76) explained an additional 4% of variance in accrual among sites. Included in the model but not clearly related to biomass accrual was days post flood and SRP (VIP scores < 0.9). The biomass sampling in June occurred just 5 days after spring flood, which might have had a scouring effect in some streams (M1 and M2, Fig. S6B) but this effect was not significant in the overall PLSR model. Finally, light was negatively related to biomass in the model, likely because it was negatively correlated with DIN.

Seasonality of GPP
Daily GPP ranged from 0 to 4.9 g O 2 m −2 d −1 , with an average among streams during the open water season of 0.8 g O 2 m −2 d −1 (Fig. 2). While there was no detectable GPP in any stream from December to January, the length of "aquatic" growing season varied considerably among sites in the network. First, the onset of measurable GPP during winter/spring transition differed by as much as 3 months among sites. For example, measurable rates were observed as early as February in one of the small birch forest streams (M16), but not until June for two tundra sites (M6 and M17). Between these end members, GPP was initiated during April at M18 and M1, and during May for M2, M10, and M9. By contrast, the end of the aquatic growing season was more constrained across the network with measurable rates of GPP until the end of November in the birch forest streams and in one tundra stream (M6), but only until mid-October for the other tundra streams.
The initiation and termination of GPP at the spring and autumn transitions were closely related to light and temperature. As light and temperature are highly correlated during these time windows, as well as temporally autocorrelated with GPP in some cases, it was not possible to separate the individual responses in GPP, yet we found higher productivity at low light in autumn compared to spring. For example, GPP at M2 and M16 ranged from 0 to 4 g O 2 m −2 d −1 in October and November when estimates of incident light were 0-150 μmol m −2 s −1 and water temperature ranged from 0.5 to 5 C. By comparison, in February and March at these same sites, GPP was only 0-0.6 g O 2 m −2 d −1 despite similar light levels, but far colder water temperatures (< 0.5 C, Fig. S8). Consistent with this, the GPP and light correlations in autumn had a much higher slope (0.27 in M16 and 0.69 in M2) compared to in spring (0 in M16 and 0.05 in M2).
For the open water season (May to October), we observed distinct temporal patterns in daily GPP among streams (Fig. 2). Four sites had no clear pattern: two tundra streams (M6 and M17) maintained notably low rates of GPP (average of 0.3 g O 2 m −2 d −1 ) throughout spring, summer, and autumn, while the two small birch forest streams (M16, M18) had comparatively high rates (average 1.5 O 2 m −2 d −1 ) throughout these time windows. The other streams in the network showed obvious periods with elevated GPP. At the main stem (M1), maximum GPP occurred early in the season; M9 (tundra) and M2 (mixed tundra/birch) had maximum GPP in autumn, whereas one tundra stream (M10) had elevated rates in both shoulder seasons.
GPP during the open water season was largely constrained by nutrient concentrations (Fig. 3, Table S1). In this case, N and P were both important predictors of daily GPP, but their relative importance differed among sites based on the background DIN concentration. In detail, daily GPP increased with DIN concentration in three streams (M10, M9, and M2, R 2 > 0.4) that all had average DIN concentrations below 30 μg L −1 . GPP for these same streams was unrelated to both temperature and light (Table S1). In comparison, GPP increased with SRP alone or in combination with light in the three streams that had DIN > 40 μg L −1 (M1, M16, and M18, R 2 = 0.3-0.6 for SRP and 0-0.5 for light). Finally, GPP was unrelated to both nutrients and light at the two streams which lacked seasonality in GPP (M6 and M17, DIN 34 and 9 μg L −1 , respectively). Considering the proportional area of each study stream in the broader network (Table 1), roughly 66% of the open water GPP varied in response to changes in N concentration alone, while 17% responded to P and only 2% to a combination of P and light.

Magnitude, patterns, and controls of annual GPP
Annual GPP varied ninefold among streams, from 47 to 411 g O 2 m −2 yr −1 (Fig. 4). Differences in annual GPP across sites were positively related to the average DIN concentration (PLSR, correlation score of 0.87, VIP score 1.57) and negatively related to winter season length (correlation score −0.91, VIP score 1.47) and the number of flood events (correlation score −0.47, VIP score 0.97, 82% of variance explained by the first two components). Slope, SRP, temperature, and drainage area were included in the analysis but were not important predictors of GPP based on VIP (scores < 0.9). Of this annual GPP, one stream (M16) had significant winter productivity of 15%. Further, a large fraction of annual GPP occurred during shoulder seasons, between 12% and 77% among sites, with M2 having the highest and M17 the lowest (Fig. 4).
To generate an estimate of seasonal GPP for the full network, we scaled up daily rates for individual streams based on their proportional areal coverage within this drainage system. The two largest streams that together make up 60% of the catchment are the tundra streams M9 and M10 (Table 1), followed by the main stem, M1 (17%). Results from this scaling exercise suggests a clear bimodal pattern of network productivity, with maximum GPP in May and then again in late August to early September, while the lowest open water GPP occurred in July (Fig. 5). When compared to estimates of terrestrial GPP in the catchment, which peak in mid-summer, this network-scale pattern of aquatic production reveals a clear temporal mismatch. At the network scale, we estimate that only 52% of annual aquatic GPP was produced during summer months/terrestrial growing season (June to August, 5B) with ca. 46% occurring during the shoulder periods (May, September, and October).

Discussion
In this first assessment of continuous, annual metabolism for Arctic or sub-Arctic streams, we show that nutrient availability is important for the seasonality of GPP and the cumulative annual productivity. In this catchment, nutrient concentrations are lowest during the terrestrial growing season, when incident light is elevated across the stream network. Consistent with our expectations, the opposing seasonality of these drivers led to a temporal coupling between nutrients and GPP and resulted in the highest productivity rates during the shoulder seasons when light availability and inorganic resources were aligned. In fact, up to 77% of the annual GPP accumulated outside the summer months. Our results further highlight the importance of ice cover duration, especially during spring, which set the length of the aquatic growing season and the potential for streams to take advantage of periods when nutrients are elevated. Finally, this study reveals a clear temporal offset in terrestrial vs. aquatic productivity in this landscape, which appears to be regulated largely by the seasonal timing of nutrient use and release on land.

Drivers of GPP seasonality
Extended periods of darkness during winter exert obvious constraints on annual patterns of stream GPP (Huryn et al. 2014). Despite this, "growing season" lengths differed among streams in the Miellajokka network by almost double, owing largely to differences in the timing of ice-off in spring. Furthermore, our observations suggest that for the streams that were open earliest, GPP could be temporarily suppressed by low temperature, although such effects are difficult to isolate without experiments (e.g., see Hood et al. 2018). Another factor which could be important in regulating the onset of GPP in early spring is carryover of viable autotrophic biomass from the previous autumn (Kendrick and Huryn 2015). We saw evidence for such carryover of biofilms in one tundra stream (M10) and bryophytes at one birch forest stream (M16), both of which incidentally had elevated rates of GPP in the spring. Regardless, our results highlight important variability in the length of the aquatic growing season within this single drainage network, and this ultimately influenced both seasonal and annual GPP patterns.
While long winters at high latitudes create uniformity in the overall GPP seasonality, patterns observed throughout the spring, summer, and autumn were notably distinct among sites. These differences during the open water season were unrelated to light but were instead explained by variation in nutrient concentrations over time. In this catchment, and consistently across the Arctic, nutrient concentrations drop from the onset of spring melt and are typically lowest in July (Holmes et al. 2000;Mann et al. 2012). In Miellajokka, the decline in DIN is greater than would be expected based solely on dilution (Figs. S4 and S5), suggesting that this seasonality is, in part, related to upregulated terrestrial nutrient demand during the growing season (Weintraub and Schimel 2005b). This pattern creates persistent nutrient limitation of algal growth during summer (Peterson et al. 1983;Pastor et al. 2020), which appeared to suppress summer GPP in all but the two most nutrient rich streams (> 90 μg DIN L −1 ) and led to positive correlations between the temporal variability in GPP and DIN concentrations at the more N poor sites (20-30 μg L −1 ). Notably, based on the relative area of different streams, DIN-driven metabolic patterns during the open water season appear to operate across ca. 60-70% of the channel network, including most of the streams draining tundra vegetation that dominates this regional landscape.
These correlations are generally consistent with the common occurrence of N as a limiting nutrient in terrestrial (Sitters et al. 2019) and aquatic (Bergström et al. 2013;Myrstener et al. 2018) ecosystems of northern Sweden. However, relatively small changes in chemistry among streams also appeared to alter the relative importance of N and P as drivers of GPP variability. In our catchment, this shift appears linked to the deep glacial sediments in the alluvial valley, which support productive birch forests and supply more DIN to streams when compared to tundra counterparts. For sites influenced by these sediments, temporal patterns of GPP were uncoupled from variability in DIN concentration and were instead correlated with either P (in the main stem) or a combination of nutrients and light (in the small birch forest streams). While it is interesting to observe differences in the role of N and P within the same drainage network, the broader message is that, for oligotrophic high latitude streams, both nutrients are important. Indeed, similar transitions between N, P, and NP co-limitation have been documented among lakes (Bergström et al. 2013) and across terrestrial ecosystems (Sundqvist et al. 2014) in this region. Thus, assumptions about the identity of a single limiting nutrient should be made with caution as this can be dynamic at small spatial and temporal scales.
Overlaying the influence of seasonal light and nutrient regimes gave rise to elevated GPP during the shoulder seasons. For sites that opened early in the spring, the time between iceoff and the spring flood represented a critical period, when primary producers could capitalize on the combination of high light and high nutrients (M1, M18, and M16). Similarly, from late summer until winter, GPP was clearly elevated at several sites (M2, M9, and M10) through increases in DIN, until light limitation eventually set in. Algal accrual and standing stock data support these seasonal patterns-and further suggest the highest accrual rates during October. Surprisingly, metabolic regimes with autumnal GPP peaks are rare in published records, but were observed in individual reaches in Methow River basin, Washington state (Mejia et al. 2019) and indicated from chamber-based estimates of GPP in the Kuparuk River, Alaska (Kendrick and Huryn 2015). While poorly studied at high latitudes, these periods of high biomass production may serve as important sources of energy to food webs prior to aquatic insect emergence in early summer, and before entering a period of strong carbon limitation during winter (Huryn and Benstead 2019).
In addition to the influences of light and nutrient supply, our results indicate that flood disturbance likely exerted both positive (Ulseth et al. 2018) and negative (Uehlinger and Naegeli 1998) effects on temporal patterns of GPP. For example, during the spring flood in 2017 and in May 2018, an initial flush of P in the absence of scouring appeared to enhance rates of GPP in the main stem (M1). By comparison, in response to a second spring flood peak in June 2018, DIN at this site dropped to 27 μg L −1 and at the same time rates of GPP decreased to summer low values (see Figs. 2 and S3). While this second flood may have directly reduced GPP (e.g., via scouring), the apparent lack of post-flood recovery, as well as the low biomass accrual sustained throughout summer were more likely a result of low nutrients than disturbance per se. For comparison, after a resetting scour event at a high DIN site (M18) during summer 2016, GPP returned to preflood level within 2 weeks (Fig. 2), similar to recovery rates observed elsewhere, including alpine streams (Uehlinger 2006). These observations suggests that rapid post flood recovery is possible in this landscape, but requires sufficient nutrient supply. Importantly, given that our data mostly encompass a single year that was relatively dry and included few large floods, our results potentially understate the role of flood disturbance in this system.

Annual and network GPP patterns
While the seasonality of GPP determines the timing of peak energy supply and biogeochemical processing, annual productivity represent the total amount of energy fixed in a system with direct consequences for consumers (Kendrick et al. 2019).
In the Miellajokka catchment, median, annual GPP was low relative to published records (Ulseth et al. 2018;Savoy et al. 2019), yet, rates at our most productive stream were very similar to the Ivishak Spring, Alaska (Huryn et al. 2014). Unlike expectations for temperate streams and rivers (e.g., Finlay 2011), variation in annual GPP was unrelated to drainage size, likely because the interactions between light and drainage size that typify forested biomes do not operate in tundra landscapes. Instead, cumulative productivity increased among streams with average N concentration and the length of the ice-free period. Separating these two drivers is difficult, however, because the most nutrient rich streams also tended to be open early in the spring. Taken together, the highest annual productivity corresponded to streams that could capitalize on summer incident light by having elevated nutrient supply and only moderate disturbance (Koenig et al. 2019;Savoy et al. 2019). In this catchment, and probably at high latitudes in general, these high GPP streams account for a small fraction of the total drainage length, yet their ecological significance for consumers is likely high (Parker and Huryn 2013).
Although we observed distinct temporal patterns among streams, together these formed a clear seasonality at the network scale, with elevated GPP in both spring and autumn, primarily influenced by the larger tundra streams which dominate this landscape. This emergent pattern also reveals a clear temporal mismatch in the phenology of terrestrial vs. aquatic GPP in this landscape (Fig. 5). Such an offset is consistent with the situation where the distribution and movement of a single resource, in this case N, constrains the productivity of two connected ecosystems. Sharing the same limiting resource on land and water creates a strong "functional coupling" in space, where the donor system (here: terrestrial plants and soils) has primacy over resource use, and thus determines its supply to the recipient (here: streams). According to this hypothesis, stream GPP is partly reliant on seasonal inefficiencies in terrestrial N retention. Similar coupling of terrestrial and aquatic GPP via material exchange has recently been suggested for lakes (Walter et al. 2020), and such relationships are potentially even stronger in running waters, given the more direct and rapid connection to catchment soils. This type of land-water interaction is distinct from situations where terrestrial constraints on aquatic productivity occur passively (e.g., via canopy shading), or not at all if hydrologic disturbance acts as the main driver of aquatic productivity. The strength and nature of such a functional connection is important because it means that any climate changes in the Arctic that upregulate terrestrial nutrient use on land (e.g., "greening"), or facilitate soil nutrient losses (e.g., permafrost thaw; Salmon et al. 2018), are likely to have immediate consequences for the productivity regimes of recipient streams and rivers.

Conclusions
Here we provide a set of mechanisms that we think likely underpin productivity regimes of streams in northernmost Fennoscandia. One likely feature of these regimes is that a large fraction of annual GPP occurs during shoulder seasons, when the timing of ice and snow cover relative to nutrient availability is crucial for the magnitude of daily stream GPP. Our results further suggest that the exact temporal patterns of GPP during the open water season will likely vary among individual streams and rivers based on the details of seasonal nutrient supply. In the context of climate change, the net effects of warming, greening, and altered hydrology on nutrient export from Arctic catchments remain unresolved (e.g., Weintraub and Schimel 2005a), but likely depend on the occurrence of thawing permafrost, which can act as a strong N source (Salmon et al. 2018). Yet, Arctic and sub-Arctic areas without significant permafrost have not received much attention even though at the latitude of our studied catchment (68 N), 70% of the global land mass lacks permafrost altogether or has only discontinuous or isolated cover (Zhang et al. 1999). Thus, Miellajokka is representative of many nutrient poor catchments in the lower Arctic and sub-Arctic that may not experience large scale nutrient increases as the climate warms.
Finally, our results call for a closer examination of the effects of nutrient supply on the productivity regimes of streams and rivers. While others have documented coupling between nutrient concentrations and daily GPP across sites or between years (e.g., (Mulholland et al. 2001;Bernot et al. 2010;Arroita et al. 2019), we show that such relationships also operate at seasonal time scales. At high latitudes, nutrient-driven patterns during the open water season are perhaps not surprising given the combination of high light conditions and strong likelihood of nutrient limitation. Yet, for any stream or river with persistently elevated light and fluctuating nutrient concentrations, resource supply could drive productivity regimes. Such regulation should only be important if N and/or P concentrations are periodically depleted to limiting levels, which may not be common in regions with widespread anthropogenic enrichment. For example, observations from larger temperate rivers show light-driven peaks in summer GPP that must be supported by sufficient nutrient supply throughout the season . Similarly, where land-use activities remove riparian canopies, the effect of nutrients supply on patterns of stream GPP could be amplified, but may also be masked by other drivers, like altered thermal and disturbance regimes (Beaulieu et al. 2013;Blaszczak et al. 2019). Regardless, we suggest that a better integration of nutrient and productivity data will increase our understanding of metabolic patterns in running waters, particularly in remote settings.