Permafrost Boundary Shift in Western Siberia May Not Modify Dissolved Nutrient Concentrations in Rivers

Identifying the landscape and climate factors that control nutrient export by rivers in high latitude regions is one of the main challenges for understanding the Arctic Ocean response to ongoing climate change. This is especially true for Western Siberian rivers, which are responsible for a significant part of freshwater and solutes delivery to the Arctic Ocean and are draining vast permafrost-affected areas most vulnerable to thaw. Forty-nine smalland medium-sized rivers (10–100,000 km2) were sampled along a 1700 km long N–S transect including both permafrost-affected and permafrost-free zones of the Western Siberian Lowland (WSL) in June and August 2015. The N, P, dissolved organic and inorganic carbon (DOC and DIC, respectively), particular organic carbon (POC), Si, Ca, K, Fe, and Mn were analyzed to assess the role of environmental parameters, such as temperature, runoff, latitude, permafrost, bogs, lake, and forest coverage on nutrient concentration. The size of the watershed had no influence on nutrient concentrations in the rivers. Bogs and lakes retained nutrients whereas forests supplied P, Si, K, Ca, DIC, and Mn to rivers. The river water temperature was negatively correlated with Si and positively correlated with Fe in permafrost-free rivers. In permafrost-bearing rivers, the decrease in T northward was coupled with significant increases in PO4, Ptot, NH4, pH, DIC, Si, Ca, and Mn. North of the permafrost boundary (61◦ N), there was no difference in nutrient concentrations among permafrost zones (isolated, sporadic, discontinuous, and continuous). The climate warming in Western Siberia may lead to a permafrost boundary shift northward. Using a substituting space for time scenario, this may decrease or maintain the current levels of N, P, Si, K, Ca, DIC, and DOC concentrations in rivers of continuous permafrost zones compared to the present state. As a result, the export flux of nutrients by the smalland medium-sized rivers of the Western Siberian subarctic to the Arctic Ocean coastal zone may remain constant, or even decrease.


Introduction
Prediction of the macro-and micro-nutrient export from the land to the ocean under various climate change scenarios is among the major scientific challenges of aquatic biogeochemistry [1][2][3][4][5].This is particularly true for high latitude regions, which exhibit large variation in soil properties, plant communities, and the distribution of aquatic systems.In this regard, Siberian watersheds draining both peatlands and mountain regions are at the forefront of scientific efforts because of their important role in element (C and nutrients) delivery to the Arctic Ocean [6][7][8][9].In particular, Western Siberian watersheds are responsible for a significant part of the freshwater and solutes delivery to the Arctic Ocean [9][10][11].These rivers flow through permafrost-affected areas which are most vulnerable to thaw [12].
The significance of high latitude continental zones in the regulation of global C and nutrient cycles stems from high storage of C org , P, and N in organic-rich soils, and from the dominance of a humid climate supporting abundant vegetation and yielding high riverine fluxes of solutes from the land to the Arctic Ocean.Climate change is likely to have a profound influence on the magnitude of these fluxes by altering the hydrological regimes of major Arctic rivers [8,13] and stimulating element release caused by widespread permafrost thaw [14].In fact, increased transport of C org , P, and N may significantly change primary productivity in riverine [15,16], estuarine [15,17,18], and ocean ecosystems [19], making predictions of climate change impact on Arctic terrestrial-aquatic ecosystems a major concern.
Due to the high importance of the Arctic Ocean and permafrost-dominated sub-arctic continental zones in the global carbon cycle and the high vulnerability of circumpolar zones to climate warming [20][21][22], numerous studies have been devoted to the biogeochemistry of organic carbon and nutrients in large rivers of the circumpolar zone [3,[23][24][25][26][27][28][29][30].The main source of information on long-term fluxes of dissolved and suspended material from Northern Eurasia to the Arctic Ocean is data from the Russian Hydrometeorological Survey (RHS) collected at key gauging stations of almost all major Russian Arctic rivers [31][32][33][34].These data comprise the daily discharge and several (typically from five to 15 per year) measurements of dissolved load.In contrast to the dissolved inorganic load (Ca, Mg, Si, Cl, SO 4 , HCO 3 ) which is adequately assessed by the RHS, for nutrients (N, P, C, Fe, Mn), a rigorous procedure of sampling, filtration, storage, and analysis is largely absent.This limits the availability of information on nutrient export from Siberian rivers solely to the data from the Great Arctic River Observatory [35][36][37] sampled several times a year at the terminal gauging stations of large Russian rivers draining to the Arctic Ocean.While this assessment allowed quantification of large Siberian river fluxes and their seasonal and spatial pattern, it provided little insight into potential environmental (landscape) controls of N and P concentration in rivers.Such factors may include the mean annual air temperature (MAAT) in the watershed, runoff, proportion of wetlands, forests and lakes, biomass productivity, permafrost extent and the depth of the active layer.Large river monitoring does not enable distinguishing between these factors and, consequently, to foresee the response of riverine nutrient fluxes to climate change.This paper aims to fill this gap by investigating small (<100-100,000 km 2 watershed area) Western Siberian rivers across a sizeable latitudinal gradient of permafrost zones.
Large continental plains, such as the Western Siberian Lowland (WSL), are especially important in nutrient transport to the Arctic Ocean, notably due to a thick peat layer [38] rich in C, N, and P. The catchments of the WSL are highly homogeneous in landscape, lithology, and topography.As a result, the permafrost, MAAT and vegetation may be the main factors controlling riverine nutrient exports to the Arctic Ocean [9,[38][39][40].Taking advantage of this unique geographical situation, a "substitution space for time" approach was successfully used in Western Siberia for inorganic river load [10], dissolved organic carbon [39,40], trace metals [41] and nutrients [9].The global picture emerging from these studies is that southern, permafrost-free regions export 3 times more solutes [10], N [9], P [8], and that wetlands exert a significant positive effect on carbon and nutrient concentrations in small rivers [8,9].However, the timing of the former sampling was mainly restricted to one-point measurement during summer baseflow conditions (August and July 1999, 2000, and 2001 [9,10,39]).While this approach is useful for distinguishing nutrient mobilization from the soil active layer and groundwater reservoirs, such sampling does not allow straightforward prediction of the change in the riverine nutrient export to the Arctic Ocean, because the majority (60% to 80%) of water and solute fluxes in high latitude zone occur during the spring flood (May-June) as shown in various Siberian regions [31,42], and the Mackenzie River [18].The purpose of this study was to improve the understanding of the magnitude and seasonality of the riverine nutrient export both in spring freshet and summer baseflow.For this, we aimed at quantifying concentrations of macro-(N, P, Si, C, K, Ca) and micro-nutrients (Fe, Mn) across a vast permafrost and landscape gradient (1700 km) with special emphasis on the permafrost-bearing zone during the two main hydrological regimes of 2015, the peak of the spring flood (June) and the summer baseflow (August).Our previous work in this region was devoted to the biogeochemistry of carbon and major and trace elements in rivers over the four main seasons of 2014 (February, June, August, and October) [40,41], whereas such nutrients as N, P, and POC have not yet been addressed before.
Our working hypothesis, in accordance with Frey et al. [9], is that northern regions should exhibit a lower concentration of nutrients both during the spring flood and summer baseflow conditions, whereas a shift of the permafrost boundary northward will enhance the nutrient export from small rivers to large rivers and subsequently to the Arctic Ocean.Our objectives were to (i) quantitatively test the role of wetland, lake, and forest coverage, MAAT and permafrost extent on nutrient concentrations in rivers draining to the largest frozen peatland territory of the world; and (ii) employing an approach substituting space for time, developed for the surface waters of Western Siberia, [9,10,[39][40][41], predict possible changes in nutrient concentrations in conditions of climate warming and progressive permafrost thaw.

Study Site and Methods
The WSL is located between the Ural Mountains in the West and the Yenisey River in the East.This vast (>1 million km 2 ), relatively pristine peatland and forest zone exhibits several physico-geographical features making it a unique, but at the same time representative, territory for testing the effect of climate, vegetation, and permafrost distribution on the hydrochemistry of nutrients in river water.The WSL is situated in the taiga forest, forest-tundra, and tundra zones.The position of natural biomes follows a decrease in MAAT from −0.5 • C in the south to −9.5 • C in the north.The peat thickness ranges from 1 to 3 m, but is not directly linked to the latitude and exhibits a uniform layer of organic-rich material [38,43].The lithology of underlining sedimentary rocks varies from silt and clay deposits with carbonate concretions in the south, to clay and sand Neocene deposits in the north.The runoff gradually increases from 100 mm in the south to 300 mm in the north [44] due to a decrease in evaporation and evapotranspiration northward, despite a sizeable decrease in precipitation from the south (500-600 mm y −1 ) to the north (360 mm y −1 ).
The 49 rivers of the WSL sampled during the spring flood (15-25 June) and summer baseflow (25 July-21 August) 2015 belong to the watersheds of the Ob, Pur, and Taz (Figure 1).The sampling strategy for the freshet was designed to sample at the peak of the spring flood.For, this we moved along the main roads from the south to the north progressively sampling the rivers, following the northward advance of the snowmelt.Note, however, that it is impossible to sample various rivers across a 1700-km length at the exact same stage of the hydrography.The physico-geographical characteristics of the catchments (Table S1 of Supplementary Information) were determined by digitalizing available soil, vegetation, lithological, and geocryological maps.Water samples were collected from the middle of the river at a depth of 0.5 m in pre-cleaned polypropylene bottles and were filtered immediately on site.Samples (50 mL) of filtered (through pre-combusted at 550 • C for 4 h acid-washed 0.45 µm Whatman GF/F glass fiber filters, (Thermo Fisher Scientific, Waltham, MA, USA) water were immediately frozen for later analysis of NH 4 + -N, NO 3 − -N, and PO 4 −3 -P.The nutrient concentrations were analyzed using an automated flow injection analyzer (FIA star 5000, FOSS, Denmark) with detection limits of 1 µg L −1 for N-NH 4 + , 0.5 µg L −1 for N-NO 3 , and 0.5 µg L −1 for P-PO 4 3− .The filters remaining after filtration were stored frozen until analysis of sestonic particulate C using an IL 550 TOC/TN analyzer with an incineration unit (1300 • C), following prior acidification to remove all inorganic C. Samples for DOC, DIC, major cations, anions, total dissolved phosphorus (P tot ), Si, Fe, and Mn were filtered through 0.45 µm acetate cellulose filters (Millipore, Sartorius Burlington, MA, USA) and analyzed as described previously [40,41].Multiple regression analysis was performed for quantifying the relationship between nutrient concentration and other physico-geographical parameters of the watershed, including MAAT, water temperature, permafrost, wetlands, lake, and forest coverage.More thorough statistical treatment of nutrient concentrations in river waters included a Principal Component Analysis (PCA) analysis using the Statistica package comprising the methods for scores and variables [45].To identify the group of nutrients that behaved in a similar way in river water, a complementary hierarchical cluster analysis (HCA) was applied (e.g., [46]).
Water 2017, 9, 985 4 of 18 C using an IL 550 TOC/TN analyzer with an incineration unit (1300 °C), following prior acidification to remove all inorganic C. Samples for DOC, DIC, major cations, anions, total dissolved phosphorus (Ptot), Si, Fe, and Mn were filtered through 0.45 μm acetate cellulose filters (Millipore, Sartorius Burlington, MA, USA) and analyzed as described previously [40,41].Multiple regression analysis was performed for quantifying the relationship between nutrient concentration and other physico-geographical parameters of the watershed, including MAAT, water temperature, permafrost, wetlands, lake, and forest coverage.More thorough statistical treatment of nutrient concentrations in river waters included a Principal Component Analysis (PCA) analysis using the Statistica package comprising the methods for scores and variables [45].To identify the group of nutrients that behaved in a similar way in river water, a complementary hierarchical cluster analysis (HCA) was applied (e.g., [46]).The Pearson correlation distance was used to trace the linkage distance, and it was considered significant at 0.1.The PCA analysis was used for the full set of sampled rivers for all seasons simultaneously, for each season individually, and separately for small rivers (<10,000 km 2 ) of the permafrost-affected zone (north of 61 • N).Here, the average latitude and area of the watershed, river water temperature, the mean annual air temperature of the watershed, runoff, permafrost coverage and percentage of bogs, forest, and lakes, the pH, and all nutrient concentrations were considered as numerical variables.In addition, a general assessment of the permafrost impact on river water nutrient concentration was performed by separating all the sampled watersheds into five categories according to the permafrost distribution in the WSL: permafrost-free (south of 61 • N), isolated (61 to 63.5 • N); sporadic (63.5 to 65 • N); discontinuous (65 to 66 • N); and continuous permafrost zones (north of 66 • N).The normality of variables was assessed using the Shapiro-Wilk test.Since the data were not normally distributed, non-parametric statistics were used.Thus, the median, first, and third quartiles were used to trace the dependence of nutrient concentration on the type of permafrost distribution.The difference in nutrient concentration between each two adjacent permafrost zones was tested using the Mann-Whitney U test for a paired data set with a significance level at 0.05.

Element Correlations, PCA, and Cluster Analysis
The pairwise correlations of component concentrations in all rivers during both seasons yielded a statistically significant relationship (R Pearson > 0.40 at p < 0.05) of PO 4 and P tot with NO 3 , DIC, Si, K, Ca, and Mn (Table 1).The N-NO 3 was not linked to any other nutrients except NH 4 .POC and DOC did not significantly (R < 0.4) correlate with any other component of the river water.Considering permafrost-free and permafrost-bearing watersheds separately, the correlation between nutrients and other solutes was mostly pronounced in rivers located south of the permafrost boundary (Table S2 of the Supplementary Information).Strong correlations of P tot and PO 4 with DIC, Si, K, and Ca remained both in permafrost-free and permafrost-bearing zones, but NH 4 did not correlate with any component (except POC) in permafrost-bearing rivers whereas it was strongly linked to PO 4 , POC, DIC, Si, Ca, and Mn in permafrost-free rivers (Table S2 of the Supplementary Information).The N-NO 3 positively correlated with PO 4 , N-NH 4 , DIC, and Ca in permafrost-free rivers, but did not exhibit any statistically significant correlations with other solutes in the permafrost-bearing zone.The PCA treatment of nutrient concentration in the WSL rivers did not identify any distinct factors controlling the overall variability in the dataset.The 13 components (pH, PO 4 , NO 3 , NH 4 , POC, DOC, DIC, Si, P tot , K, Ca, Fe, and Mn) and eight physico-geographical parameters of the watershed (latitude, area, runoff, MAAT, percentage of bogs, forest, and lakes and the permafrost coverage) were tested during all seasons together and separately for June and August.However, this did not yield more than three factors with a variability explanation of ≤6% for the first two factors and ≤3% for the third one.Considering solely small rivers (<10,000 km 2 watershed area) did not improve the explanation capacity of the PCA.
The hierarchical cluster analysis did not yield any significant link between most components at the distances <0.1, regardless of seasons (Figure S1A of the Supplementary Information).Only in August, K, Ca, DIC, and NO 3 − demonstrated some links at ~0.1.During this season, small rivers located north of 61 • N, at the permafrost boundary, demonstrated significant links between PO 4 , K, and NO 3 at ~0.1 (Figure S1B of Supplement 1).

Latitudinal Pattern of Nutrient Concentration
The pairwise correlations of component concentrations and landscape parameters including the latitude of the watershed are presented in Table S3 of the Supplement.Generally, considering all rivers and all seasons simultaneously (Table 2), there was no correlation of latitude with all solutes except DOC (R Pearson = −0.48).In permafrost-free rivers, the latitude negatively correlated with P tot and PO 4 , whereas in permafrost-bearing rivers, the latitude negatively correlated only with DOC and K (Table S3).In June, there was no effect of latitude on PO 4 , P tot , NO 3 , NH 4 , POC, and Si concentration (R 2 ≤ 0.1, p > 0.05).In August, there was no effect of latitude on PO 4 , P tot , NO 3 , NH 4 , POC, DOC, and Si whereas DIC, K, and Ca exhibited increasing concentration northward (not shown).
To synthetize the observed latitudinal trends within a context of permafrost distribution, five main permafrost zones were selected: absent (south of 61 • N), isolated (61 to 63.5 • N), sporadic (63.5 to 65 • N), discontinuous (65 to 66 • N) and continuous (north of 66 • N).A box plot (median, first, and third quartile) of nutrient concentrations during spring and summer in each permafrost zone (Figure 2) demonstrated that the majority of macro-and micro-nutrients did not show statistically significant dependence on permafrost coverage.Only DIC, K and Ca had significantly (p < 0.05) higher concentrations in the permafrost-free zone compared to permafrost-bearing zones.The difference in nutrient concentration between each two adjacent permafrost zones was tested using the Mann-Whitney U test for a paired dataset (Table S4 of the Supplementary Information).The main result is that, both in June and August, there was no difference in PO 4 , P tot , NO 3 , NH 4 , and POC concentration in four pairs of zones: continuous/discontinuous; discontinuous/sporadic; sporadic/isolated; and all permafrost-bearing sites/no-permafrost zone.Only pH and concentrations of DIC, K, and Ca were lower in the permafrost-bearing zone compared to rivers located south of the permafrost boundary.The pH, DOC, DIC, Ca, and Mn were also different between continuous and discontinuous permafrost zones but the concentrations of DOC, DIC, Ca, and Mn were actually higher in the continuous permafrost zone compared to the discontinuous one.The spring season yielded the minimal differences between permafrost zones: among all studied nutrients, only K was higher in permafrost-free zones compared to permafrost-bearing ones.In August, there was an enrichment in DOC, DIC, NO 3 , Si, K, and Ca in permafrost-free zones compared to permafrost-bearing zones, whereas pH, Ca, and Mn concentrations were higher in the continuous permafrost zone compared to the discontinuous one.

Impact of Landscape and Physico-Geographical Parameters on Individual Components
The impact of the several most important landscape parameters (permafrost, forest, bog coverage, and river water temperature) on nutrient concentration was tested using Pearson correlation for all rivers (Table 2) and separately for permafrost-bearing and permafrost-free zones (Table S3 of the Supplementary Information).P tot and PO 4 demonstrated a positive correlation with forest coverage, visible only in permafrost-bearing watersheds, whereas NO 3 positively correlated with forest coverage only in permafrost-free watersheds.The concentration of POC did not show any statistically significant relationship with all the environmental factors in both permafrost-free and permafrost-bearing zones.
In contrast to N and P, concentrations of other nutrients in rivers were to a various degree related to landscape parameters.The degree of permafrost coverage negatively correlated with DOC, Si, K, Ca, Mn, and Fe (Table S3, Figure S2 of the Supplementary Information).The increase in runoff raised the concentration of PO 4 , DIC, K, and Ca and decreased the concentrations of NH 4 and Mn in permafrost-bearing watersheds.In contrast, the permafrost-free watersheds did not demonstrate any effect of runoff on nutrient concentration (Table S3).The watershed area positively correlated with K (R = 0.57) in permafrost-free rivers, but had a minor impact on other nutrients (Table 2).Indeed, the Spearman rank correlations between the watershed area and pH, DOC, DIC, PO 4 , P tot , NO 3 , NH 4 , POC, Si, K, Ca, Fe, and Mn were not significant for all rivers (Table 2) or in permafrost-free rivers only (Table S3).Considering both seasons and all rivers simultaneously, only permafrost-bearing watersheds yielded a statistically significant Pearson correlation between S watershed and pH, K, and Ca concentrations.Moreover, the impact of the watershed area (S watershed ) on nutrient concentration was tested using non-parametric Spearman rang correlation coefficient (Table S5 of the Supplementary Information).Only small rivers (S watershed < 100 m 2 ) showed a positive correlation of S watershed with pH, DIC, PO 4 , P tot , POC, Ca, and Mn in June and with pH, DIC, PO 4 , P tot , K, Ca, Mn, and Fe in August, whereas the watershed area of larger rivers did not control the nutrient concentrations.
River water temperature (T river ) negatively correlated with DIC, Si, and Mn for all the rivers of the dataset (Table 2), and only with Si in permafrost-free rivers.In permafrost-bearing rivers, the increase in river temperature was coupled with a significant (R > 0.22 at p < 0.05) decrease in pH, DIC, Si, NH 4 , PO 4 , P tot , Ca, and Mn (Table S3 of the Supplementary Information).The watershed air temperature (T watershed ) positively correlated with DOC for the whole dataset (Table 2), with PO 4 and P tot in permafrost-free watersheds and with DOC and K in permafrost-bearing rivers (Table S3).Other nutrients did not exhibit any statistically significant link with watershed air temperature.
The bogs decreased pH, DIC, Si, and Mn in all WSL rivers (Table 2).In permafrost-free watersheds, bogs decreased NO 3 , DIC, Si, K, and Ca concentrations, although the correlations were weak (Table S3).In permafrost-bearing rivers, bogs also decreased pH, PO 4 , P tot , DIC, Si, K, Ca, and Mn, as illustrated in Figure 3.The lakes were negatively correlated with PO 4 , P tot , DIC, Si, K, Ca, and Mn concentrations in permafrost-bearing rivers (Figure 4).The forest exerted a positive effect on pH and the concentration of NO 3 , DIC, and Ca in permafrost-free watersheds, and on pH, PO 4 , P tot , DIC, Si, K, Ca, and Mn in permafrost-bearing ones (Table S3).The effect of the forest on NO 3 , NH 4 , PO 4 , Si, Ca, and Fe concentration during June and August across the full latitudinal profile of the WSL is illustrated in Figure 5.Note that the percentage of lakes + bogs and forest was negatively correlated.

Discussion
The obtained results will first be compared with previous data on nutrients in the WSL rivers, then the impact of environmental parameters on the variability of nutrient concentration across the latitudinal and permafrost gradient will be analyzed and, finally, the possible impact of the permafrost boundary shift northward on nutrient concentration in WSL rivers will be discussed.
We believe that the period of spring flood sampling chosen in this study is representative for all the rivers of the transect.First, a particularity of the Western Siberian Plain is that the spring freshet period is much longer than in the majority of other subarctic settings.Very flat lowlands, the abundance of bogs, and thick peat coverage lead to a broad peak of water flow, which lasts more

Discussion
The obtained results will first be compared with previous data on nutrients in the WSL rivers, then the impact of environmental parameters on the variability of nutrient concentration across the latitudinal and permafrost gradient will be analyzed and, finally, the possible impact of the permafrost boundary shift northward on nutrient concentration in WSL rivers will be discussed.
We believe that the period of spring flood sampling chosen in this study is representative for all the rivers of the transect.First, a particularity of the Western Siberian Plain is that the spring freshet period is much longer than in the majority of other subarctic settings.Very flat lowlands, the abundance of bogs, and thick peat coverage lead to a broad peak of water flow, which lasts more than one month.Analyses of stable water isotopes in soil waters, lakes, and rivers in the Western Siberian Lowland on a large spatio-temporal scale along a 1700 km transect [47] demonstrated that (i) a small proportion of the water which contributed to rivers flow during the spring was released from melting snow; and (ii) the primary runoff generation mechanism in the WSL is the displacement of water already stored in the lakes, wetlands, and soils to rivers as shallow subsurface flow and suprapermafrost flow.For this reason, the exact timing of the freshet and the position of the sampling point at the river hydrograph are of secondary order importance for understanding the nutrient pathways in the riverine systems.Moreover, the variation in hydrology is believed to play a limited role in DOC (and, presumably, nutrients) variability and export from the WSL watersheds [40].Nevertheless, part of the variability in the dataset can be contributed by the rapidly changing hydrograph (especially in northern rivers), the varying freshwater sources along the south-to-north gradient (especially in spring), and the fact that it is impossible to sample 40 rivers at the exact same stage of the hydrography.
The concentrations of Si, Ca, DIC, DOC, K, Mn, and Fe in the WSL rivers sampled in 2015 in this study are in fair agreement with those measured in 2013 and 2014 within the same latitudinal profile [40,41].The average P tot , NO 3 , and NH 4 concentrations in the WSL rivers in summer (104 µg P L −1 ; 40 µg N-NO 3 L −1 , and 59 µg N-NH 4 L −1 ) are higher (Mann-Whitney U test at p < 0.0003) compared to previous studies in the region by Frey et al. [9] (48 µg P L −1 , 35 µg N-NO 3 L −1 , and 7.6 µg N-NH 4 L −1 , respectively) for a similar subset of rivers in terms of latitudinal and permafrost zone position.In agreement with Frey et al. [9], our results show no statistically significant difference in NO 3 between permafrost-influenced and permafrost-free watersheds although the difference in NH 4 concentration is clearly pronounced.The NO 3 and NH 4 concentrations of the WSL rivers are in general agreement with other small rivers of the permafrost zone in Alaska [26].McClelland et al. [37] reported 1.2-2.4ppm POC in the Ob River, in fair agreement with the values 0.5-2.5 ppm in this study for WSL rivers including the Ob River middle course.
In the permafrost-free rivers sampled in August, the concentration of NO 3 , NH 4 , DIC, Ca, and K decreased northward, which is in general agreement with former observations on N [9] and inorganic nutrients [8,10].A systematic northward decrease in Ca, DIC, and nutrients (N, Si, K, and Mn) in permafrost-free watersheds is interpreted as being due to the influence of shallow underground and subsurface waters that interact with the mineral substrate and feed the river during the open water period.The chemical composition of this main feeding source systematically changes across the latitudinal profile.The presence of soluble minerals in background substrates is also traced by Sr and Mo concentration in WSL rivers [41].In the southern part of the WSL, carbonate concretions and limestone-rich clays that underlay peat and forest soils [40,48] may enrich the river waters in P, Ca, DIC, and Mn.In accord with this general lithological context of the WSL, a dramatic difference (by a factor of ca.[5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] in DIC and Ca concentrations between permafrost-free and permafrost-bearing zones was recently interpreted as due to the presence of highly-soluble carbonate concretions in clay-silt bedrocks which underlie forest and peat soils in the south of the WSL [40].In addition to this lithological impact, in summer, high concentrations of NO 3 , NH 4 , Si and K in some rivers of this region (see Figure S2 of the Supplementary Information) may be due to nutrient leaching from the soil litter of the taiga forest, abundant in the southern watersheds.
A striking result of this study was the lack of effect of key environmental "climatic" parameters such as the latitude and degree of the permafrost coverage on most nutrient concentrations in rivers situated north of the permafrost boundary.Moreover, a number of components (PO 4 , P tot , NH 4 , Si, DIC, Ca, and Mn) demonstrated a significant (at p < 0.05) increase in concentration with a decrease in water temperature in permafrost-affected watersheds, and DOC, Si, K, Ca, Mn, and Fe was negatively correlated with permafrost coverage of the watershed (Figure S2 of the Supplementary Information).In the permafrost zone, the mineral substrate does not interact with suprapermafrost water that feeds the river [49].Unlike in southern rivers, shallow subsurface and active layer waters interact exclusively with peat layers in the northern watersheds.It is thus possible that the high homogeneity of the peat layer in the northern part of the WSL creates a low contrasting pattern of nutrient concentration in rivers.Namely, peat leaching may be the main limiting factor of this nutrient release in permafrost-affected watersheds.In contrast, in permafrost free-rivers, only Si negatively correlated with water temperature, which can be explained by an enhanced biological uptake of Si by plankton, periphyton and macrophytes in the river channel with the increase in water temperature.
The lakes retained nutrients such as P, Si, K, and Mn (Figure 4), as is well known across the globe [50].The retention of major and micronutrient export from watersheds by lakes and bogs (wetlands) is also well known [51] and generally consistent with our observations on PO 4 , P tot , DIC, Si, Ca, K, and Mn across the WSL territory.Positive correlations of PO 4 and P tot with DIC, Si, K, Ca, and Mn (see Table 1 and Table S2 in Section 3.1) suggest a common mineral source of these elements in rivers, since P, unlike N, is mostly derived from chemical rock weathering [52,53].Indeed, the highly positive impact of the forest coverage of the watershed on nutrient (PO 4 , Si, P tot , K, Mn), CO 2 , DIC, and Ca concentrations in the permafrost-bearing zone may be due to feeding the rivers by shallow subsurface waters.These waters are present only in forested zones, typically along the river valley, where the permafrost is absent or its thickness is particularly low [40,47,54,55].In such settings, the groundwaters that feed the river interact with mineral horizons and can be enriched in "mineral" nutrients such as P, Si, and K, whereas Mn may be mobilized in reduced form from Mn hydroxide nodules, abundant at the base of the sand horizon [48,56].The effect of forest presence is more pronounced in the permafrost-bearing regions because the surrounding frozen peat can release a lesser amount of nutrients than the deep mineral soil horizons in permafrost-free regions.Alternatively, the positive effect of the forest coverage on nutrient concentration in permafrost-bearing watersheds may be due to recycling of tree litter that is capable of releasing a sizeable amount of these biogenic elements [57,58].Unlike in essentially forested permafrost-free regions, the presence of even a small fraction of forest along the river valley in the generally tundra and palsa context of the northern part of the WSL (e.g., [59,60]) produces a visible impact on nutrient concentration in the river water.The forest litter leaching by surface waters may be responsible for the positive relationship of NO 3 and NH 4 with the percentage of forest coverage in permafrost-free watersheds (Figure 5), given that the N, unlike P, is not sensitive to mineral weathering [61].Note that in Finnish catchments, subjected to excess atmospheric N loading, the concentration of NH 4 was higher in the catchments with a high peatland proportion compared to those with a low peatland proportion [62,63].
The lack of a watershed size (S area ) impact on nutrient concentration in rivers (Table 2) suggests highly homogeneous landscape settings across a wide range of watershed size, from 100 to 100,000 km 2 .Only in rivers <100 km 2 watershed area, a positive effect of S area on pH, DIC, P, K, Ca, Mn, and Fe (Table S5 of the Supplementary Information) may be due to non-homogeneity of their watersheds, i.e., an increase in forest versus bog proportion with the increase in S area .In very small watersheds, the feeding of streams occurs via leaching of surface peat.As such, the appearance of shallow groundwater draining some mineral layers under the permafrost-free forest zone may produce a marked effect on nutrient concentrations.Presumably, the S area > 100 km 2 presents a certain threshold after which the non-homogeneity of the watershed disappears and the effect of groundwater becomes independent from the size of the catchment.
As can be seen from Figure 2, the spring-time DOC was much higher than summer-time DOC in the permafrost-free boreal zone, which is consistent with general understanding of DOC behavior in subarctic catchments.The difference between seasons in the permafrost-bearing zone is quite low, which may be due to strong dilution by melted snow and the lack of forest in permafrost-bearing zones.Another reason for this behavior is the peatland context of all the sampled watersheds.According to the results of our previous studies [40,49], in the permafrost zone, the DOC (and, presumably, nutrient) export is limited by water travel time through organic topsoil, litter, and ground vegetation (moss, lichen) rather than through the peat and mineral layer.In that case, it is only the thickness of the unfrozen peat and the local permafrost coverage that control the DOC export from the soil to the river.As a result, the DOC and nutrient concentration in the streams will be weakly dependent on the watershed size and seasons.
It is known that autochthonous processes within the river channel are capable of retaining nutrients [64] and metals [65].Indeed, with a rather low annual runoff of the WSL rivers, sizeable retention of dissolved N, P, Si, Fe, Mn, and Al by oxyhydroxides and Si by coastal grass and diatoms in the river may occur.However, given that the size of the river (and, thus, the water residence time in the channel) have an insignificant effect on the nutrients considered in this study (Table S5 of the Supplementary Information), we argue that there is a negligible impact of in-stream retention on nutrient transport in WSL rivers.
Overall, highly-contrasting patterns of latitudinal impact on nutrients in permafrost-free and permafrost-bearing zones were observed.The main novel finding is that there was only weak or no latitudinal dependence on nutrient concentration in the permafrost zone and the difference between adjacent permafrost zones was not significant not only during the summer baseflow, but also during the spring flood.The lack of a permafrost coverage effect on N, P, and POC concentrations in rivers located north of the permafrost boundary may require revising the predictions of climate change impact on WSL river biogeochemistry and nutrient transport to the Kara Sea.The concept of "substituting space for time" in the WSL hypothesized in earlier works [9,10,39] suggests that there will be a shift of the permafrost boundary to the north, leading to a progressive replacement of the type of permafrost distribution (continuous → discontinuous; discontinuous → sporadic; sporadic → isolated, etc.) rather than global transformation of permafrost-bearing territory to permafrost-free territory.The data presented in Figure 2 and statistical tests of the difference between adjacent permafrost zones (Table S4 of the Supplementary Information) demonstrate that, during the open water period (spring + summer), there will be no significant (at p < 0.05) change of P tot , PO 4 , NO 3 , NH 4 , DOC, POC, K, Si, DIC, Fe, and Mn concentrations in the case of a permafrost boundary shift and the change of the permafrost type distribution.
A tentative explanation for this is that permafrost thaw facilitates export during the summer in all permafrost zones, but in southern rivers a greater pool of nutrients is available for export compared to the northern ones.However, this is largely compensated by more rapid flushing and a shorter travel time through soils and rivers, as well as lower microbial activity in rivers of the continuous permafrost zone.Thus, contrasting gradients in supply vs. in stream removal may cancel out the net effect of temperature and permafrost on nutrient concentrations in the river water.

Conclusions
The impact of watershed parameters (latitude, temperature, bog, lake, and permafrost coverage) on main nutrient (N, P) dynamics in WSL rivers during spring flood and summer baseflow was found to be quite low.Considering all watersheds during both periods, only the proportion of forest coverage was positively correlated to the concentrations of PO 4 , P tot , NO 3 , Si, K, and Mn.A two-to 10-fold decrease in DIC, Ca, and K concentrations from the permafrost-free zone to the beginning of permafrost appearance was presumably linked to the change of the underlying rock lithology and the flow paths from underground and shallow subsurface flow to surface and active layer flow.
Overall, the nutrients exhibited a highly contrasting pattern of latitudinal impact: weak or no latitudinal effect in the permafrost zone, and a northward decrease in the permafrost-free zone.In rivers located north of the permafrost boundary, significant negative relationships between water temperature and nutrient concentrations were observed for P, Si, DIC, Ca, and Mn, whereas an increase in permafrost coverage decreased the concentrations of Si, K, Ca, Mn, and Fe, but did not impact the concentrations of PO 4 , P tot , NO 3 , NH 4 , POC, and DIC.It is possible that, despite lower pools of nutrients in northern soils compared to permafrost-affected regions, more rapid flushing and a shorter travel time through soils and rivers, as well as lower microbial activity in rivers of the continuous permafrost zone, provide the concentrations of nutrients similar to those in discontinuous and sporadic permafrost regions.
Applying the concept of substituting a space for time scenario in WSL rivers, we foresee that the permafrost boundary shift may not modify most nutrient concentrations greater than 30-50% due to the change of continuous to discontinuous permafrost on a long-term scale in the most northern rivers.Therefore, we do not anticipate any dramatic impact of air temperature rise and progressive permafrost thaw on nutrient export from soils to rivers and, subsequently, to the Arctic Ocean.

Figure 1 .
Figure 1.Sampling sites and physico-geographical context of WSL territory investigated in this work.

Figure 1 .
Figure 1.Sampling sites and physico-geographical context of WSL territory investigated in this work.

Figure 2 .
Figure 2. Box plot of first and third quartiles (25% and 75%) of nutrient concentrations in five permafrost zones during August (white rectangles) and June (grey rectangles).The horizontal line represents the median and the whiskers represent the minimum and maximum value.The outliers are not shown.

Figure 2 .
Figure 2. Box plot of first and third quartiles (25% and 75%) of nutrient concentrations in five permafrost zones during August (white rectangles) and June (grey rectangles).The horizontal line represents the median and the whiskers represent the minimum and maximum value.The outliers are not shown.

Figure 3 .
Figure 3. Correlations (R is the Pearson coefficient) between K, Ca, DOC, Fe, and bog proportion on the watershed in permafrost-free (open circles) and permafrost-bearing (closed diamonds) rivers.The two sampling periods, June and August, are combined together.The solid line is for the main trend and the dashed lines are significant at p < 0.05.

Figure 4 .
Figure 4. Si, DIC, K, and Mn relationships (R is the Pearson coefficient) with lake proportion on the watershed in permafrost-free (open circles) and permafrost-bearing (closed diamonds) rivers sampled in June and August.The two sampling periods, June and August, are combined together.The solid line is for the main trend and the dashed lines are significant at p < 0.05.

Figure 3 .
Figure 3. Correlations (R is the Pearson coefficient) between K, Ca, DOC, Fe, and bog proportion on the watershed in permafrost-free (open circles) and permafrost-bearing (closed diamonds) rivers.The two sampling periods, June and August, are combined together.The solid line is for the main trend and the dashed lines are significant at p < 0.05.

Figure 3 .
Figure 3. Correlations (R is the Pearson coefficient) between K, Ca, DOC, Fe, and bog proportion on the watershed in permafrost-free (open circles) and permafrost-bearing (closed diamonds) rivers.The two sampling periods, June and August, are combined together.The solid line is for the main trend and the dashed lines are significant at p < 0.05.

Figure 4 .
Figure 4. Si, DIC, K, and Mn relationships (R is the Pearson coefficient) with lake proportion on the watershed in permafrost-free (open circles) and permafrost-bearing (closed diamonds) rivers sampled in June and August.The two sampling periods, June and August, are combined together.The solid line is for the main trend and the dashed lines are significant at p < 0.05.

Figure 4 .
Figure 4. Si, DIC, K, and Mn relationships (R is the Pearson coefficient) with lake proportion on the watershed in permafrost-free (open circles) and permafrost-bearing (closed diamonds) rivers sampled in June and August.The two sampling periods, June and August, are combined together.The solid line is for the main trend and the dashed lines are significant at p < 0.05.

Figure 5 .
Figure 5. N-NO3, N-NH4, P-PO4, Si, Ca, and Fe relationships with forest proportion on the watershed during June (open circles) and August (closed diamonds) in all rivers of the WSL sampled in this work.The R value is the Pearson correlation coefficient for the full dataset.The solid line is for the main trend and the dashed lines are significant at p < 0.05.

13 Figure 5 .
Figure 5. N-NO 3 , N-NH 4 , P-PO 4 , Si, Ca, and Fe relationships with forest proportion on the watershed during June (open circles) and August (closed diamonds) in all rivers of the WSL sampled in this work.The R value is the Pearson correlation coefficient for the full dataset.The solid line is for the main trend and the dashed lines are significant at p < 0.05.

Table 1 .
Correlation matrix of dissolved components (D.C.) and nutrient concentrations in all rivers sampled during June and August.Marked (bold) Pearson correlations R > 0.40 are significant at p < 0.05 (N = 74).

Table 2 .
Correlation matrix of watershed physico-geographical parameters and dissolved component (D.C.) concentrations in all rivers sampled in June and August.Marked (bold) Pearson correlations R > 0.40 are significant at p < 0.05 (N = 74).S is the area of the watershed; T river is water temperature at the day of sampling; T watershed is mean annual temperature of the watershed.Lat and Permaf represent latitude ( • N) and permafrost coverage of the watershed, %.The bogs, forest, and lakes represent the percent coverage in the watershed.