Browning affects pelagic productivity in northern lakes by surface water warming and carbon fertilization

Abstract Global change impacts important environmental drivers for pelagic gross primary production (GPP) in northern lakes, such as temperature, light, nutrient, and inorganic carbon availability. Separate and/or synergistic impacts of these environmental drivers on pelagic GPP remain largely unresolved. Here, we assess key drivers of pelagic GPP by combining detailed depth profiles of summer pelagic GPP with environmental and climatic data across 45 small and shallow lakes across northern Sweden (20 boreal, 6 subarctic, and 19 arctic lakes). We found that across lakes summer pelagic GPP was strongest associated with lake water temperatures, lake carbon dioxide (CO2) concentrations impacted by lake water pH, and further moderated by dissolved organic carbon (DOC) concentrations influencing light and nutrient conditions. We further used this dataset to assess the extent of additional DOC‐induced warming of epilimnia (here named internal warming), which was especially pronounced in shallow lakes (decreasing 0.96°C for every decreasing m in average lake depth) and increased with higher concentrations of DOC. Additionally, the total pools and relative proportion of dissolved inorganic carbon and DOC, further influenced pelagic GPP with drivers differing slightly among the boreal, subarctic and Arctic biomes. Our study provides novel insights in that global change affects pelagic GPP in northern lakes not only by modifying the organic carbon cycle and light and nutrient conditions, but also through modifications of inorganic carbon supply and temperature. Considering the large‐scale impacts and similarities of global warming, browning and recovery from acidification of lakes at higher latitudes throughout the northern hemisphere, these changes are likely to operate on a global scale.


| INTRODUC TI ON
Global environmental changes such as warming, recovery from acidification and forestry have altered the biogeochemistry in northern lakes Skjelkvåle et al., 2001) and have thus impacted important environmental drivers for pelagic gross primary production (GPP). As phytoplankton are crucial providers of energy, minerals and biochemical compounds for higher consumers, understanding global change impacts on pelagic GPP is highly relevant for aquatic food webs in general (Müller-Navarra, 2008;Peltomaa et al., 2017;Sterner & Hessen, 1994). At higher latitudes, lakes are particularly common (Lehner & Döll, 2004;Verpoorter et al., 2014) and surface air temperature anomalies related to climate change are the greatest (the arctic amplification) (Cohen et al., 2014;Hansen et al., 2010;Serreze et al., 2009), emphasizing the need of understanding global change impacts on northern lakes in particular. Increased air temperatures induce a direct warming of surface waters (O'Reilly et al., 2015;Schneider et al., 2009), but also indirectly increase precipitation Hudson et al., 2003;Lind & Kjellström, 2008). Warming and increased precipitation further promote enhanced catchment vegetation cover (i.e., greening), which together with forestry and recovery from acidification induce enhanced loadings of terrestrial dissolved organic material (DOM) to northern lakes Finstad et al., 2016;Kritzberg, 2017). Important components of DOM related to lake biogeochemistry are nutrients and colored dissolved organic carbon (DOC). However, neither the separate nor the synergistic impacts of large-scale changes of warming and altered lake biogeochemistry on pelagic GPP at higher latitudes are resolved.
Pelagic GPP is commonly measured at discrete depths (here GPP z rates) resulting in differently shaped depth profiles depending on light availability and the maximum GPP rate (here: GPP z,max , per m 3 ) occurs where optimal growth conditions are present (Wetzel & Likens, 1991;Figure 1a). The GPP z rates can be upscaled to a lakeaverage (GPP lake-average , per m 2 ) by integrating rates over the water column and dividing them by the lake surface area. Key environmental drivers therefore likely differ between GPP z,max rates that depend on local conditions and GPP lake-average , which reflects the response to integrated environmental conditions.
The pelagic GPP z rates can be constrained by dissolved inorganic carbon (DIC), nutrients, light (energy), and temperature (Graham et al., n.d.;Wetzel & Likens, 1991). Carbon dioxide (CO 2 ) is the most bioavailable DIC source used in photosynthesis by phytoplankton, and pH regulates the amount of CO 2 relative to (bi)carbonates in lake water DIC (Huisman et al., 2018;Wetzel & Likens, 1991). Hence, lake pH is important for GPP via its effect on the relative amount of CO 2 .
The degree of CO 2 supersaturation of surface waters often increases with DOM concentrations in lakes (Del Giorgio et al., 1999;Larsen et al., 2011;Sobek et al., 2003). Although several mechanisms can cause CO 2 supersaturation in lakes, terrestrial DOM is important as it is correlated with lower pH in lakes and its mineralization generates CO 2 (Lazzarino et al., 2009;Nydahl et al., 2020;Stets et al., 2017).
Several empirical and modeling studies show that the tradeoff between light and nutrient availability promotes a unimodal distribution of pelagic GPP lake-average with increased lake DOC (Rivera Vasconcelos et al., 2018;Solomon et al., 2015), where the peak in GPP lake-average is determined by the DOC:nutrient stoichiometry, whereas the location of the GPP peak along the DOC axis is determined by the light climate Kelly et al., 2018). Nonetheless, the unimodal relationships vary and are not always observed, suggesting that other factors influence this relationship and regulate pelagic GPP lake-average (Kelly et al., 2018). For example, co-limitation by CO 2 and nutrient on pelagic GPP has been demonstrated in experimental studies (Low-Décarie et al., 2011, and field studies in northern oligotrophic lakes (Brown et al., 2019;Hamdan et al., 2018;Jansson et al., 2012). In addition, results from controlled experimental pond ecosystems further suggest that warming alone may additionally amplify pelagic GPP lake-average at all levels of lake water DOC concentrations (Hamdan et al., 2018;Figure 1b). Increased colored DOC also results in epilimnion warming, especially in small and shallow lakes, likely through intensified water column stratification (Bartosiewicz et al., 2016(Bartosiewicz et al., , 2019Houser, 2006;Pilla et al., 2018Pilla et al., , 2020. Since a warmer climate promotes higher lake DOC concentrations via indirect effects occurring in the lakes catchment (Laudon et al., 2012;Tetzlaff et al., 2013), increasing air temperatures can amplify the epilimnion warming via increasing DOC, that is, additional DOC-induced warming of the lake epilimnia (here named: internal warming). However, intrinsic effects of temperature on pelagic GPP lake-average may be hard to disentangle from colored DOC and nutrients (Bergström et al., 2013;Faithfull et al., 2011;Klug, 2005), and CO 2 (Jansson et al., 2012;Jonsson et al., 2001) as they are tightly correlated, making temperature redundant in models that include DOC, nutrients, and CO 2 .
Moreover, these recognized patterns and relationships between environmental drivers of pelagic GPP (rates and lake-average) may not be consistent across biomes and over seasons due to variable DOC:nutrient stoichiometry, coloring of the DOC, differences in climate, catchment properties, air temperatures, and light conditions Isles et al., 2021;Kelly et al., 2018;Seekell et al., 2015). Overall, several studies have assessed how DOC impacts light and nutrient availability and quality, lake water DIC concentrations, and temperature (summarized in Figure 1b), but the relative contribution of DIC, CO 2 , temperature, and DOC:DIC stoichiometry as additional drivers for pelagic GPP lake-average has been far less assessed in empirical studies.
Here, we investigate how global change by its impact on key environmental drivers affects summer pelagic GPP in northern lakes over a large spatial scale. For this reason, we collected data on pelagic GPP depth profiles and GPP lake-averages , together with physico-chemical lake parameters in summer (June-August) from 45 lakes in northern Sweden. The lakes were spread over three different biomes (20 boreal, 6 subarctic, and 19 Arctic lakes; Figure 1c), covering a colored DOC gradient from 1.0 to 19.5 mg L −1 . Besides abiotic in situ conditions and lake bathymetry, we also investigate the additional effects of air temperature on the DOC-induced warming of the lake epilimnia (i.e., internal warming), and how the inorganic carbon sources (CO 2 vs. DIC) are a function of lake water pH. We assess how key environmental drivers related to global change affect summer GPP (rates and lake-averages) both across northern Sweden and per biome, and how DOC:nutrient (DOC:TN; DOC:TP) and DOC:DIC stoichiometry, and lake water temperature influence GPP lake-average . We hypothesize that: (1) Internal warming of lakes increases along the DOC gradient, (2) GPP z,rates relates to temperature, DOC, nutrients DOC:nutrient stoichiometry and CO 2 , and GPP lake-average relates to similar drivers but mostly to lake bathymetry and (3) drivers are biome specific.

| Study area, sampling and data compilation
We compiled data of pelagic GPP, together with water chemistry and bathymetry, from 45 small and shallow lakes (lake surface area between 0.6 and 9.4 ha, and depth between 3.7 and 15.8 m) F I G U R E 1 (a) Conceptual figure showing depth profiles of GPP z rates in a stratified brown lake high in DOC (black), in a lake intermediate in DOC (gray) and in a clear lake low in DOC (blue). In brown lakes GPP is confined to a shallower epilimnion, but rates are high, and in clear lakes GPP is spread deeper in the water column, but with consistently lower GPP z rates. (b) Conceptual figure of pelagic GPP lake-average distribution (green) with increasing DOC concentrations, initially limited by nutrients, and at higher DOC concentrations by light inhibition. At higher DOC concentrations increased temperatures (red arrow) can increase maximum GPP z rates and could thus to some extent counteract the negative impact of reduced light on pelagic GPP. The height of the peak in GPP is defined by DOC:nutrient stoichiometry, and the location of the GPP peak by the coloring of DOC. (c) Sampling locations, with the lakes from the Arctic biome in the northernmost outlined area, subarctic in the western outlined area, and boreal in the southernmost outlined area. Map lines delineate study areas and do not necessarily depict accepted national boundaries.  Figure S1). Data and sampling methods of pelagic GPP and water chemistry for 35 of the lakes are described in detail in earlier publications (Ask et al., 2009;Deininger et al., 2017Deininger et al., , 2019Karlsson et al., 2001), and the remaining 10 lakes (four in the boreal and six lakes in the subarctic biome, original to this study) are sampled using similar techniques. The data can be found at https://doi. org/10.5061/dryad.t1g1j wt5k.

| Pelagic GPP
Pelagic GPP was measured in situ at the surface and at subsequent 1 m depth intervals, with additional measurements at 0.25 and 0.5 m, where the deepest measurement depended on the lake depth and water turbidity (sensu Ask et al., 2009;Deininger et al., 2017Deininger et al., , 2019Karlsson et al., 2001). Measurements were done by incubating (single or replicates, see raw data) transparent glass bottles filled with water from the sampling depth, with additional incubations in dark bottles at the most shallow and deepest measurements, for about 4 h midday using 14 C isotopic tracer as described by Schindler et al. (1972).
We used raw isotopic activity values from all datasets, and recalculated GPP values similarly for consistency among the datasets. We further estimated the consumption of the DIC pool in the bottles to rule out that DIC became limiting during the 4 h incubation by calculating the DIC consumption from the GPP z max in the lake over the incubation time of 4 h in relation to the DIC pool in the incubation bottle. This showed that in all cases the consumption of DIC in the bottle was less than 5% of the DIC pool (median value 0.3%), except for Struptjärn (which we excluded from the analyses) where the consumption was clearly higher and ca 62% of the pool in the bottle was consumed. Calculated GPP values were extrapolated to daily GPP using the ratio of incident photosynthetically active radiation (PAR) during incubation in relation to daily PAR. GPP z,max rates represent the maximum rate measured over the water column occurring at one specific depth in each lake (i.e., the peak in the vertical distribution; Figure 1a). GPP rates measured at discrete depths (GPP z rates) were upscaled to a single average GPP estimate per lake (GPP lake-average ) by integrating the rates over the water column and dividing them by the lake surface area (Table S1). We used volume-weighted pelagic GPP z rates when upscaling to GPP lake-average .

| Water chemistry
Since most of the pelagic GPP takes place in the epilimnion, and nutrients within the epilimnion are well mixed (all lakes stratified except for nine Arctic lakes: see raw data), we only used values from the epilimnion (i.e., 1 m) when relating water chemistry concentrations to pelagic GPP z,max rates and GPP lake-average (see earlier publications for details on sampling procedures). Samples for pH, DOC, DIC, total nitrogen (TN), and phosphorous (TP) were thus taken at 1 m depth (epilimnion) or were taken from composite water samples (lakes from Ask et al. (2009)). In short, DOC was filtered through a 0.45 μm filter (Sarstedt Filtropur), acidified using HCl to an end concentration of 12 mM, and stored in a refrigerator before analyzed. TN and TP (unfiltered) samples were kept frozen until analysis. Dissolved inorganic carbon (DIC) samples were taken by injecting 4 ml (air free) lake water into a tightly sealed 18 ml glass vial (pre-flushed with N 2 ) con-  (Isles, 2020). pH was measured immediately after sampling in the laboratory, and CO 2 concentrations in the lake water were calculated from DIC, pH, and temperature, following guidelines from the Water Quality Analysis

| Light, temperature, and bathymetry
PAR and temperature in all lakes were measured using a handheld probe at every meter throughout the water column at the deepest part of the lake, with additional measurements at 0.25 and 0.5 m.
Light attenuation coefficients (Kd) were calculated as the absolute slope of natural logarithmically transformed PAR against depth.
From the Kd we calculated the percent of incoming light at the depth where the pelagic GPP z,max was located (%light) relative to the surface (100%), and for each lake we calculated the euphotic depth (the depth where 1% of surface light remains; z euph ). Together with incoming daily PAR (defined as PAR) we also calculated the daily PAR at depth for the GPP z,max rates (PAR depth ). Daily PAR was collected from stations located next to the lake, or acquired from open databases (Laudon et al., 2021;SMHI, website). Pelagic GPP z,max rates are related to their depth specific temperature (T depth ), whereas pelagic GPP lake-average estimates are related to temperature at a fixed epilimnion depth of 0.2 m (T water ). Average air temperatures (T air ) are obtained using monthly air temperature averages 1 month before sampling from weather stations located within a maximum range of 60 km from the sampling sites (data extracted from the Swedish Meteorological and Hydrological Institute [https://www.smhi.se/ klima t/klima tet-da-och-nu]), including a temperature decrease of 0.57°C per 100 m elevation difference between station and sampling site (sensu Karlsson et al., 2005 and references therein).
Internal warming of the epilimnion is here defined as T water -T air . We calculated the lake-average depth (z avg ) and volumes (as a whole, or in different sections) and lake surface areas (Area) from detailed bathymetry (Table S1).

| Statistical analyses
Differences in water chemistry variables among the biomes are tested in ANOVAs with estimated marginal means compared per biome (Searle et al., 1980). To relate internal warming to DOC concentrations, we first selected which curve type (linear, logarithmic, or exponential) best explained (highest R 2 ) T air and T water versus DOC, and then identified the main drivers of internal warming using a multiple linear regression (MLR) with forward selection that included the variables DOC, Kd, Area, z max , z avg , z euph , and lake altitude. The significance of internal warming against DOC was tested with regression analyses.
We investigated which drivers explain GPP (GPP z,max rates and GPP lake-average ) best throughout northern Sweden (i.e., including all lakes) in an MLR with forward selection. We tested data for underlying assumptions of parametric tests, but data were not corrected for collinearity. To overcome collinearity issues, unequal sample sizes per biome, we further investigated the spread of explanatory variables and GPP in a PLS using the plsr package (Mevik & Wehrens, 2007) in R, again both including data from all lakes and in addition per biome. A PLS is a comparable method to the more wellknown principal component analysis (PCA), but specifically suitable for datasets with high predictor values compared to observations, like our dataset. VIP scores are a measure of how substantially the variable adds to the model, and the loading describes the correlation intensity between the variable and the predictor. We considered a VIP <0.9 as minimum value for a variable to substantially add to the model (Mehmood et al., 2012). Included variables for both GPP z,max rates and GPP lake-average are presented in Table 1, and we included similar variables in the MLR and the PLS.
Lake 13 (no T water data) and lake 15 (ice on lake) were removed from the dataset assessing internal warming. We removed these two lakes, together with lake 4 (DIC values two standard deviations above averages), and Struptjärn (GPP z,max and GPP lake-average two standard deviations above averages because of an invasive G. semen bloom: see Deininger et al., 2019) from the dataset including GPP. GPP z,max rates and GPP lake-averages were log-transformed to meet conditions for performing parametric tests. The quality of the multiple linear regressions models was tested using Akaike's information criterion (AIC). We verified that the warming pattern is statistically independent of the sampling year and timeframe ( Figure S1). Statistical analyses were conducted in IBM spss statistics v. 26, and in R. We considered an effect statistically significant at p < .05 (two-tailed for the ANOVAs and MLR).

| Water chemistry and internal warming
Our study lakes cover a wide range of DOC, DIC, and nutrient concentrations (Table S2) DOC concentrations increased with T air (Figure 2a; p < .000).
The relationships between temperature (in water and air) and

| Summer pelagic GPP
Pelagic GPP z,max rates including lakes from all biomes increased with and were best explained by CO 2 (56.2%) and water temperature at depth (additional 19.6%) in the MLR (Figure 3; Table 1: Regression II).
These results were confirmed in the PLS regression including lakes from all biomes, as CO 2 and DOC explained the variance in pelagic GPP z,max rates best in both the first component (39.6%) and the second component (22.0%) (Figure 4a; Table S4). Variables that were included in the model that did not substantially contribute to any of the two components (VIP <0.9; Table S4) are shown in gray in the PLS figures (Figure 4). However, different drivers of pelagic GPP z,max rates were selected per biome (Figure 4b-d). In the boreal lakes, CO 2 and DIC were selected as best drivers explaining the variance in pelagic GPP z,max rates for both the first component ( Table S4b). Although variables were quite spread out in the PLS plot, indicating slight correlation among variables, T air was tightly coupled to T water in the first component, as were lake altitude and lake surface area but to a lesser extent. The PLS plot revealed that within the first component, Kd was clustered with the drivers T air , T water , and nutrients (DOC, TN, TP, and DOC:nutrient stoichiometry to a lesser extent), whereas lake area and altitude had a negative impact and thus were negatively correlated with GPP lake-average (Figure 4e; Table 2). Moreover, both euphotic depth, lake surface area, and lake altitude had a high VIP score (1.3, 1.3, and 1.1, respectively), indicating that these variables contributed substantially to the model, and thus explained part of the variance in pelagic GPP lake-average . The PLS plots also visualize that nutrients and DOC:nutrient stoichiometry did not substantially add to the model best explaining GPP lake-average , and that they were clustering together against the other variables.
Also, for pelagic GPP lake-average , different drivers were selected per biome (Figure 4f GPP lake-average increased with CO 2 in the boreal and Arctic biome but had a decreasing trend in the subarctic (Figure 5a; negative loadings in PLS: Table S4b). Moreover, the fraction of DIC that is CO 2 was very variable among our lakes and differed significantly (p < .05) between the Arctic and boreal and was on average lowest in the Arctic (22.7%), followed by the boreal (52.8%) and then by the subarctic biome (77.1%; Figure 5b).

| Nutrient stoichiometry and pelagic GPP
GPP z,max rates increased with DOC (Figure 5c), whereas the pelagic GPP lake-average tended to be unimodally related with DOC, and most lakes had DOC concentrations below the GPP peak occurring at DOC concentrations around 9 mg·L −1 (Figure 5d). We also investigated if DOC:nutrient stoichiometry (here DOC:TN, DOC:TP) and DOC:DIC modified the GPP lake-average relationship with DOC and categorized these into different levels similar to the procedure in Kelly et al. (2018) for DOC:TP ( Figure 6). As GPP in the Arctic was clearly lower compared to the other biomes, we marked GPP in the arctic biome as a category of its own. To eliminate the effect of variable coloring of DOC on a lake's light climate, we used Kd instead of DOC when graphically referring to the DOC gradient among lakes (see Figure S2 for relationship between DOC and Kd). For small and shallow lakes the light extinction coefficient (Kd) represents a good proxy for light availability in the mixed layer (Jones, 1998;Kelly et al., 2018). Generally, GPP lake-average followed a unimodal distribution with Kd, where higher DOC:TP and DOC:DIC ratios were related to lower GPP lake-average (Figure 6a,b; for the relationship with DOC:TN see Figure S3). Our dataset was, however, not large enough to perform statistical analyses on the different stoichiometry categories. Interestingly, the relationship between GPP lake-average and the different DOC:nutrient and DOC:DIC categories were affected by T water , and GPP lake-average showed a clear peak at T water temperatures around 20°C (Figure 6c,d; Figure S3b). Again, our dataset was too small here to make statistical inferences.

| DISCUSS ION
In our study, we found empirical support for elevated warming of browner surface waters, especially in shallow lakes, and for that CO 2 availability directed by a wide range of lake water pH is an important driver of pelagic summer GPP in northern lakes. Our results point out that global warming and increases in colored DOC may affect summer pelagic GPP, not only through changes in light and nutrient availability, but also via effects on water temperature and CO 2 availability.

| Internal warming of lake water
The internal warming (T water -T air ) of surface waters was magnified at low to intermediate DOC concentrations, but was dampened at DOC concentrations >14.0 mg·L −1 (Figure 2) and decreased with average lake depth (z avg ). More specifically, internal warming decreased with 0.96°C for every increasing m in average lake depth, with an internal warming being >5.0°C in lakes with z avg <4 m (Table 1: Regression I; Figure 2c). This suggests that increasing air temperatures combined with ongoing browning of lake waters will promote the most pronounced internal warming in small, shallow, and relatively clear water lakes, that is, the most abundant lake type in northern Sweden (see in . Our results support recent research advances pointing out that warming of surface waters and consequent cooling of bottom waters of stratified lakes is related to altitude, latitude, and browning, and is stronger in lakes with smaller surface area (Bartosiewicz et al., 2019;Ficker et al., 2017;Pilla et al., 2018Pilla et al., , 2020. However, our results suggest that instead of lake area, average lake depth (z avg ) predicts internal warming best, and that even among shallow lakes (up to 16 m depth in our dataset) lake bathymetry impacts internal warming. In the context of ongoing global change, internal warming of lake surface waters is therefore likely to continue to increase not only in Sweden, but also on a global scale (IPCC, 2021;O'Reilly et al., 2015). At northern latitudes where small and shallow lakes are globally most abundant (Verpoorter et al., 2014),

F I G U R E 3
Graphical and mathematical visualization of GPP z,max rates and GPP lake-average with their units, and included variables and outcomes of both the multiple linear regression (MLR) and partial least squares regression (PLS). Variables marked with a are log-transformed and abbreviations are as follows: Area, lake surface area (hectare); DOC, dissolved organic carbon (mg·L −1 ); DIC, dissolved inorganic carbon (mg·L −1 ); PAR, daily incoming PAR at surface (W·m −2 ); PAR depth , daily incoming PAR at depth (W·m −2 ); T air , average air temperature of the previous month (°C); T depth , temperature at depth (°C); T water , temperature at 0.2 m depth (°C); TN, total nitrogen (μg·L −1 ); TP, total phosphorus (μg·L −1 ); z avg , average lake depth (m); z euph , euphotic depth (m); z max , maximum lake depth (m). and warming and lake browning are especially pronounced Pagano et al., 2014;Roulet & Moore, 2006;Solomon et al., 2015), internal warming of lake surface waters is likely to accelerate and thus affect pelagic GPP.

| Pelagic GPP in lakes over the northern Swedish landscape
Interestingly, we found that summer pelagic GPP z,max rates foremost were associated with CO 2 concentrations (Figures 3-5), and secondly with water temperature (T depth ; Figure 3). Moreover, in the PLS regression DOC was selected as second driver, and GPP z,max rates increased with DOC on the first component ( Table S5; Isles et al., 2021;Kelly et al., 2018;Seekell et al., 2015). Hence, when and how nutrient and light conditions relative to CO 2 availability and temperature conditions control pelagic GPP rates is a delicate balance and includes interactive effects. In our dataset, summer GPP z,max rates were mostly related to lake CO 2 concentrations, indicating carbon fertilization.
Regardless, although measured in the season and depth with optimal light and temperature conditions, the GPP z,max rates were still impacted by both the dampening (light) and enhancing (nutrients, temperature) effects related with DOC.
In addition, the integrated effect of various drivers on the production (GPP lake-average ) showed the strongest relation with CO 2 and secondly with T air (Figures 3, 4e and 5a). An additional part of the variance in GPP lake-average decreased with DOC:DIC, and increased with DOC:TP and TP, supporting that the peak in GPP lake-average was moderated both by DOC:nutrient and DOC:DIC stoichiometry (Table 1 and Figure 1a). Hence, or results suggest co-fertilization by CO 2 and nutrients on pelagic GPP, similar to results of previous field

F I G U R E 4 Partial least squares (PLS) regression biplots of environmental variables and (a-d) pelagic GPP z,max rates and (e-h)
GPP lake-average . The biplots include (a, e) all lakes, and lakes from the (b, f) boreal, (c, g) subarctic, and (d, h) Arctic biome. Variables adding substantially to the PLS model with a VIP >0.9 are plotted in black, and variables that do not substantially add to the model with a VIP <0.9 are plotted in gray.

All biomes
Boreal SubarcƟc ArcƟc  (Jansson et al., 2012;Vogt et al., 2017). There are several possible explanations and mechanisms behind independent co-fertilization of CO 2 and nutrients on pelagic GPP-for example, enzymes active in photosynthesis may not be fully saturated at CO 2 levels close to equilibrium and/or that the carbon concentrating mechanism of phytoplankton cells is downregulated under nutrient limited conditions (see Jansson et al., 2012 and references therein). Independent co-fertilization of CO 2 and nutrients on GPP is also shown in several experimental settings (Hamdan et al., 2018;Low-Décarie et al., 2011. Our observation, showing that the peak of the unimodal distribution of pelagic GPP lake-average with increasing light attenuation (kd) is influenced by both DOC:nutrient and DOC:DIC stoichiometry ( Figure 6a,b) supports the idea that significant changes in pelagic GPP lake-average due to increased DIC (CO 2 ) and nutrient availability only occur when there is sufficient light throughout the water column to promote pelagic GPP. Higher DOC:TP ratios are generally related to lower GPP lake-average Isles et al., 2021;Kelly et al., 2018), which is opposite to our results in the MLR. In our dataset, DOC:DIC explained variation in the peak of GPP lake-average better than DOC:TP (or the DOC:TN). Possibly, GPP lake-average was less related to DOC:nutrient compared to DOC:DIC due to different phytoplankton nutrient limitation regimes with more P to NP co-limited conditions in the subarctic and Arctic lakes, and strict N-limited conditions in the boreal lakes (Bergström et al., 2013Isles et al., 2020 and Figure 4). Interestingly, the imply that lake inorganic carbon (CO 2 , DIC) availability and temperature are additional drivers of summer pelagic GPP (rates and lakeaverages), besides light and nutrients conditions governed by DOC, in lakes across the northern Swedish landscape.

| Pelagic GPP per biome
We found slightly different drivers for GPP among biomes, likely as an effect of biome differences in climate, catchment properties, and atmospheric N deposition (Elser et al., 2009;Lewis, 2011). For the boreal biome, lakes spanned a wide gradient in DOC concentrations, and were overall browner, warmer, and richer in nutrients, with more variable pelagic GPP (rates and lake-averages) that was mainly related to CO 2 and DOC:DIC (Figure 4b,f; Table S4). The subarctic lakes had less variable lake water temperature, DOC, and nutrient concentrations, and intermediate pelagic GPP being mostly related to temperature (Figure 4c,g; Table S4). The Arctic lakes were generally of low DOC, clear, cold, nutrient poor, and some not stratified, and generally low pelagic GPP z,max rates mostly associated with DIC and TP (Figure 4d,h; Table S4). The observed low pelagic GPP in the Arctic lakes (especially at high altitudes) is suggested to be an effect of low water temperatures leading to low growth rates even during conditions of nutrient enrichment (Bergström et al., 2013). Yet, as temperatures and nutrients are all low, relatively small changes in environmental conditions are likely to promote shifts from one limiting factor to another. Furthermore, the differences in environmental drivers of pelagic GPP that we identified for the northern Swedish landscape, and for the different biomes, emphasize the importance of including different biomes when upscaling to understand climate change effects on larger scales.

| Inorganic carbon
Altogether, our results underline the importance of DIC and CO 2 for pelagic GPP in addition to nutrients, light and temperature condi-  fects on pelagic GPP are studied empirically (Hessen et al., 2017;Jansson et al., 2012;Vogt et al., 2017). The fraction of CO 2 in the DIC pool is a function of pH, slightly moderated by temperature (see Section 2.3), and is relevant for GPP since the (bi)carbonate part of DIC is generally not as favorable for phytoplankton growth (with the exception for cyanobacteria) as CO 2 (Huisman et al., 2018;Wetzel, 2001). We included both DIC and CO 2 in the models for explaining summer pelagic GPP z,max rates and GPP lake-averages, which were both strongly related to CO 2 (explaining 56% and 53% of the variance in GPP in the MLR, respectively), and GPP lake-averages additionally decreased with increasing DOC:DIC stoichiometry. In the Arctic, however, where lakes were more alkaline (higher pH) and the fraction of CO 2 was lower (Figure 5b), DIC (and not CO 2 ) was selected as a major driver (Table S1). Hence, since the contribution of CO 2 to total DIC can be highly variable when assessing lakes in different landscapes with different catchment properties (Figure 5b), DIC might be a poor proxy for CO 2 availability.
CO 2 supersaturation is ubiquitous throughout lakes at higher latitudes caused by high input of CO 2 and organic material from land (Åberg et al., 2010;Sobek et al., 2003). However, natural variability of pH across landscapes related to catchment characteristics and in situ metabolism and biological engineering (Huisman et al., 2018;Paerl et al., 2016;Verspagen et al., 2014), but also anthropogenic influences such as increasing atmospheric CO 2 concentrations and ongoing recovery from acidification (Garmo et al., 2014;IPCC, 2021;Isles et al., 2018;Skjelkvåle et al., 2001) are likely to impact the amounts and form of inorganic carbon available for photosynthesis. Enhanced inorganic carbon availability in the form of CO 2 will promote summer pelagic GPP according to our results, but could also potentially impact phytoplankton community composition, favoring species that lack the ability to use (bi)carbonate as an inorganic carbon source (e.g., chrysophytes; Bhatti & Colman, 2005;Maberly et al., 2009), relative to species that have the ability to do so (e.g., cyanobacteria; Huisman et al., 2018;Verspagen et al., 2014).
Enhanced pH accompanied with lake warming, may thus further promote cyanobacteria over other species of phytoplankton (Huisman et al., 2018;Verspagen et al., 2014).

| Drivers of pelagic productivity in lakes across the northern landscape
Our results, together with others, can be applied to identify potential trajectories for pelagic GPP in northern latitude lakes following global change. Adding to previous studies recognizing the importance of DOC, and DOC:nutrient stoichiometry, impacting light and nutrient availability Kelly et al., 2018;Isles et al., 2021;Rivera Vasconcelos et al., 2018), our study further emphasizes the role of inorganic carbon, DOC:DIC stoichiometry, and temperature for pelagic GPP in summer (summarized in Figure 7). Our results imply that for northern Sweden global warming combined with browning likely enhances summer pelagic GPP via effects on water temperature, nutrients, and CO 2 . Considering the large-scale impacts and similarities of global warming and browning on lakes at higher latitudes throughout the northern hemisphere (Cohen et al., 2014;Creed et al., 2018;Solomon et al., 2015), these recognized changes are likely to operate on a global scale.
However, our results also emphasize that caution should be made when upscaling and using space-for-time substitutions. We use summer data from different years representing only a snapshot in time, where the lake environmental conditions and subsequent responses in GPP across northern Swedish landscapes have been influenced by atmospheric acid deposition and differences in catchment vegetation cover (see Isles et al., 2018Isles et al., , 2020. Yet, when analyzing air temperature data for all years within the 20-year sampling period, apparent differences in climate (air temperature) are consistently greater among than within biomes ( Figure S1).
The sampling occasions were within the natural variation in each biome, indicating that the temperature effect observed in our results is unrelated to individual sampling years. Thus, here a spacefor-time approach is valid for assessing key environmental drivers of summer pelagic GPP across the northern Swedish landscape.
However, more research is required, and consideration needs to F I G U R E 7 Conceptual overview of how our results (text and arrows in red) add to current knowledge (text and arrows in black/gray) regarding drivers of summer pelagic GPP (green) with increasing DOC. Consensus is that GPP is initially limited by nutrients, and at higher DOC by light inhibition. The height of the peak in GPP is defined by the DOC:nutrient ratio, and the location of the GPP peak by the coloring of DOC (light climate). Our results point out additional major drivers of summer pelagic GPP, which are CO 2 concentrations (determined by lake pH and the total DIC pool) increasing with DOC that together with DOC:DIC stoichiometry further enhance GPP. Increased water temperatures with increasing DOC can promote GPP rates even more and thus partly counteract the negative impact of reduced light conditions caused by DOC on pelagic GPP.