Trade-offs Between Light and Nutrient Availability Across Gradients of Dissolved Organic Carbon Lead to Spatially and Temporally Variable Responses of Lake Phytoplankton Biomass to Browning

Northern lakes are experiencing widespread increases in dissolved organic carbon (DOC) that are likely to lead to changes in pelagic phytoplankton biomass. Pelagic phytoplankton biomass responds to trade-offs between light and nutrient availability. However, the influence of DOC light absorbing properties and carbon–nutrient stoichiometry on phytoplankton biomass across seasonal or spatial gradients has not been assessed. Here, we analyzed data from almost 5000 lakes to examine how the carbon–phytoplankton biomass relationship is influenced by seasonal changes in light availability, DOC light absorbing properties (carbon-specific visual absorbance, SVA420), and DOC–nutrient [total nitrogen (TN) and total phosphorus (TP)] stoichiometry, using TOC as a proxy for DOC. We found evidence for trade-offs between light and nutrient availability in the relationship between DOC and phytoplankton biomass [chlorophyll (chl)-a], with the shape of the relationship varying with season. A clear unimodal relationship was found only in the fall, particularly in the subsets of lakes with the highest TOC:TP. Observed trends of increasing TOC:TP and decreasing TOC:TN suggest that the effects of future browning will be contingent on future changes in carbon–nutrient stoichiometry. If browning continues, phytoplankton biomass will likely increase in most northern lakes, with increases of up to 76% for a 1.7 mg L−1 increase in DOC expected in subarctic regions, where DOC, SVA420, DOC:TN, and DOC:TP are all low. In boreal regions with higher DOC and higher SVA420, and thus lower light availability, lakes may experience only moderate increases or even decreases in phytoplankton biomass with future browning.

carbon-phytoplankton biomass relationship is influenced by seasonal changes in light availability, DOC light absorbing properties (carbon-specific visual absorbance, SVA 420 ), and DOC-nutrient [total nitrogen (TN) and total phosphorus (TP)] stoichiometry, using TOC as a proxy for DOC. We found evidence for trade-offs between light and nutrient availability in the relationship between DOC and phytoplankton biomass [chlorophyll (chl)-a], with the shape of the relationship varying with season. A clear unimodal relationship was found only in the fall, particularly in the subsets of lakes with the highest TOC:TP. Observed trends of increasing TOC:TP and decreasing TOC:TN suggest that the effects of future browning will be contingent on future changes in carbon-nutrient stoichiometry. If browning continues, phytoplankton biomass will likely increase in most northern lakes, with increases of up to 76% for a 1.7 mg L -1 increase in DOC expected in subarctic regions, where DOC, SVA 420 , DOC:TN, and DOC:TP are all low. In boreal regions with higher DOC and higher SVA 420 ,

INTRODUCTION
The highest concentration of lakes in terms of both total number and area is located in northern latitudes (Verpoorter and others 2014). Many northern lakes have experienced recent increases in dissolved organic carbon [DOC; hereafter used as an analog for dissolved organic matter (DOM)] driven by a combination of changing climate, recovery from acidification, and changes in landcover and land-use activities (Monteith and others 2007;Weyhenmeyer and Karlsson 2009;de Wit and others 2016a;Kritzberg 2017). The relationship between DOC and phytoplankton biomass has been shown to exhibit trade-offs between light and nutrient availability at different concentrations of DOC (Ask and others 2009;Seekell and others 2015a;Bergströ m and Karlsson 2019). This tradeoff exists because increases in colored DOC (sometimes referred to as lake browning) may have a negative influence on phytoplankton biomass because light absorption by DOC suppresses photosynthesis and constrains growth (Carpenter and others 1998;Jones 1998;Klug 2002;Deininger and others 2017), but a positive fertilization effect because phosphorus (P) and nitrogen (N) are cotransported with DOC and these organically complexed nutrient pools frequently constitute a large proportion of total nutrient inputs to lakes (Kortelainen and others 2006;Hessen and others 2009;Kissman and others 2016).
The trade-off between light and nutrient availability may, as shown in empirical (Hanson and others 2003;Ask and others 2009) and theoretical (Kelly and others 2018;Vasconcelos and others 2019) studies, result in a unimodal response in pelagic primary production (PP) to increases in colored DOC, with peaks in PP observed at intermediate DOC concentrations of 10 to 12 mg L -1 . However, it is not yet clear under what conditions a unimodal response is likely to occur, and what implications this has for the future trajectory of PP with future browning. Several studies investigating the relationship have focused on the DOC-PP relationship, but the relationship of DOC to phytoplankton biomass has received less attention. PP and phytoplankton biomass are closely related across large gradients and aggregated timescales (Smith 1979), but the relationship is far from perfect, in part because phytoplankton biomass is affected by both growth and death processes, and high biomass may be sustained with relatively low rates of PP if grazing and other loss rates are low. Like PP, phytoplankton biomass has also been suggested to have a unimodal relationship to DOC others 2018, Bergströ m andKarlsson 2019). Theoretical models suggest that the height of the peaks in DOC-phytoplankton biomass or DOC-PP curves is related to the nutrient stoichiometry associated with DOC (Kelly and others 2018), and to the availability of multiple limiting nutrients in the case of co-limited systems (Bergströ m and Karlsson 2019). Because the hypothesized mechanism for the unimodal response involves a trade-off between light and nutrient availability, the degree to which DOC absorbs light in the photosynthetically active spectrum should also influence the position of the unimodal curve, affecting the DOC concentration at which peak phytoplankton biomass occurs (Kelly and others 2018). These theoretical models have been tested in several experimental studies, but have not often been examined using field datasets.
Here, we used data from almost 5000 lakes sampled by Swedish national monitoring programs to examine if phytoplankton biomass as indicated by chlorophyll-a (chl-a) concentration varied spatially with increasing total organic carbon (TOC) (our proxy for DOC). The lakes in this dataset span a large range, from clear and nutrient poor montane and subarctic lakes in the north and west, to brown and nutrient poor boreal lakes in central Sweden, to browner but more nutrient rich boreal lakes in the south of Sweden. The range of TOC represented in these datasets cover the range of DOC reported in global lake surveys (Sobek and others 2007). Our study builds on earlier studies in similar lakes assessing the effects of atmospheric N deposition on lake nutrient limitation in Sweden (Isles and others 2018), the effect of browning on the relative importance of N versus P limitation (Isles and others 2020), the response in pelagic primary production and phytoplankton biomass in whole-lake fertilization experiments across a DOC gradient (Bergströ m and Karlsson 2019), and the effects of browning on the stoichiometry and nutritional content of lake seston for upper trophic levels (Bergströ m and others 2020). We analyzed trends in TOC color and carbon-to-nutrient stoichiometry using both long-term and spatially resolved monitoring datasets. We tested the following hypotheses: (1) the TOC concentration at which peak chl-a occurs varies seasonally as a function of changing light availability (due to lower surface irradiance and deeper mixed layers later in the year), with the peak occurring at a higher TOC in the spring (high light) and at lower TOC in the autumn (low light); (2) the TOC concentration at which the peak chl-a occurs varies with carbonspecific light-absorption (SVA 420 ), with higher SVA 420 resulting in maximum chl-a at lower DOC concentrations; and (3) the peak chl-a varies with TOC-nutrient stoichiometry, with the height of the peak chl-a increasing with lower TOC:N or TOC:P stoichiometry. After addressing these hypotheses, we assessed long-term trends in variables that may influence the TOC-chl-a relationship. We then developed a regression model to predict future changes in epilimnetic chl-a concentrations in different parts of Sweden with an additional increase in TOC of 1.7 mg L -1 (the median change from 1996-2018) to assess the response of chl-a in regions with different lake physical and biogeochemical properties to continued lake browning.

Lake Data
The long-term monitoring (LTM) lakes were used to examine the presence of a unimodal relationship between TOC and chl-a, and to explore factors that influence the DOC versus chl-a relationship, including effects of seasonal changes in light availability (hypothesis 1), light-absorbing properties of organic matter (TOC-specific visible light absorbance at 420 nm, SVA 420 ; hypothesis 2), and TOC-nutrient stoichiometry (TOC:TP and TOC:TN; hypothesis 3). The LTM dataset includes 110 rela-tively pristine lakes sampled multiple times each year collected under the ''Trend Lakes'' and ''Reference Lakes'' programs (Supplemental Figure 1; h ttp://miljodata.slu.se/mvm/). Chl-a, total organic carbon (TOC), absorbance at 420 nm (A 420 , measured on a 5 cm path length and converted to units of m -1 )), TN and TP data from 1996 to 2018 were used in this study (exception: TN data were only available beginning in 2007). SVA 420 was calculated as A 420 /TOC for concurrent observations. Lakes with at least 10 years of data were used in the study. Lake data were divided into three seasons for analyses: April-June for the spring season; July-August for the summer season; and September-October for the autumn season. A broader time window was used for the spring season to account for the different timing of the spring bloom and the corresponding availability of monitoring data in the different regions in Sweden.
The Cyclically Low-Intensity Monitored Lakes survey (CLIML) dataset was used to assess spatial patterns in SVA 420 and carbon-nutrient stoichiometry and to predict likely future changes in chl-a with continued browning. The CLIML dataset contains 4800 lakes (4768 contained data on all variables and were used in this study) sampled once in a six-year cycle (800 per year) at autumn overturn. Data on TOC, SVA 420 , TN, and TP from the CLIML survey conducted between 2007 and 2012 were used. Chl-a data were not available in the CLIML dataset.
Samples collected from lakes in both the LTM and CLIML monitoring programs were processed in the same laboratories using similar methods (Fö lster and others 2014). Two key assumptions were made with the LTM and CLIML data. First, we assumed TOC was equivalent to DOC; the LTM and CLIML programs collect consistent data for TOC but not DOC. Previous Scandinavian studies showed that more than 90% of TOC is in the DOC pool (Kortelainen and others 2006;Weyhenmeyer and Jeppesen 2009); therefore, we assume that TOC was equivalent to DOC in our interpretation of the results, although this assumption was likely invalid in a few lakes where high nutrients resulted in high phytoplankton biomass (and correspondingly high particulate carbon). Second, we assumed chl-a was an effective proxy for phytoplankton biomass. Although chl-a content per unit biomass can vary depending on species, light conditions, and nutrient limitation (Nicholls and Dillon 1978;Kruskopf and Flynn 2006), chl-a has been shown to be a good proxy for biomass in Swedish lakes (Bergströ m and Jansson 2006), particularly across large gradients as considered in this study.

Estimation of Areal Chl-a
Theoretical models describing the pelagic DOC-PP relationship often consider areal estimates (epilimnetic depth-integrated) of PP, and studies focusing on PP often measure it in areal units based on, for example, surface oxygen and/or C 14 measurements (Hanson and others 2003;Ask and others 2009). Because lake mixed layer depth decreases with increasing dissolved organic matter (Fee and others 1996), areal units are more likely to exhibit a unimodal relationship (hypothesis 1), and estimates of areal chl-a are more appropriate for comparison with much of the literature on the DOC-PP relationship. However, most lake monitoring programs do not collect areal chl-a concentrations; rather they collect discrete epilimnetic chla concentrations which are not corrected for thermocline depth (hereafter referred to as volumetric concentrations). Therefore, for the LTM lakes, in addition to using measured volumetric concentrations, we also estimated areal chl-a concentrations (mg m -2 ) by multiplying modeled thermocline depth (m) by measured chl-a concentration (mg m -3 ). Thermocline depths were estimated using a regression model built from data from 27 sampled lakes distributed from southern to northern Sweden that were collected in 2016 (see Isles and others 2020 for a description of the sampled lakes, and Supplemental Figure 1 for their locations). The regression model incorporated lake area and DOC concentration as predictors (as in Fee and others 1996; Kelly and others 2018) and had a R 2 = 0.61 and p = 0.001: log 10 ðthermoclineDepthÞ ¼ À0:02 DOC þ 0:27 log 10 ðlakeAreaÞ þ 0:37 Thermocline depth was log-transformed to avoid predictions of negative thermocline depths in small lakes with high DOC. DOC was estimated as 0.9 9 TOC, as in Weyhenmeyer and Karlsson (2009). Our model for estimating thermocline depth assumed that all lakes were stratified, and our estimate of areal biomass assumed that hypolimnetic phytoplankton biomass was negligible. Thermal profiles to calculate thermocline depth were only available during the summer months, so this analysis was restricted to the summer season.

Statistical Analyses
To test for the presence of a unimodal relationship between TOC and chl-a, a linear function (y = a + bx), a quadratic function (y = a + bx + cx 2 ), and a Gaussian function (y ¼ ae À xÀb ð Þ 2 = 2c 2 ð Þ ) were fit to the median chl-a concentration using nonlinear quantile regression (tau = 0.5) with the function nlrq in the R package ''quantreg'' (Koenker 2013). Maximum likelihood of each curve was estimated assuming a log-normal distribution of residuals around the predicted value, and used to estimate Akaike's information criterion (AIC) for each model using the formula AIC = 2 k -2(ln L), where k is the number of free parameters and L is the loglikelihood. The model with the lowest AIC was chosen as the best fit. This method may introduce some non-independence among consecutive samples taken in the same lake; however, temporal semivariogram analysis (Supplemental Figure 9) indicated that most variance in the chl-a data in each lake could be accounted for within 1-year lags, so the effects of temporal autocorrelation on our regression results were assumed to be minimal. This method of curve fitting and model selection was employed for the entire dataset, and for the dataset subset to spring, summer, and fall seasons, to test the hypothesis that the TOC concentration at which peak chl-a occurs varies seasonally as a function of changing light availability (hypothesis 1). In the fall, data were screened to include only samples with TOC no greater than 35 mg L -1 preserving 99.7% of the original data; quadratic regression fits gave negative results beyond this threshold that could not be fit to a lognormal error distribution.
Additionally, median chl-a, as well as interquartile ranges between the 25th and 75th percentile of chl-a, was calculated for the LTM data for the April-October period using equally spaced bins of TOC and plotted, to visually assess the occurrence of a unimodal relationship. Bins with less than 10 samples were excluded from plotting. This method accounts for uneven numbers of samples at low vs. high concentrations of TOC and does not make any assumptions regarding the true functional relationship between TOC and chl-a.
To test the hypotheses that the TOC concentrations at which peak chl-a occurs vary with SVA 420 (hypothesis 2) and TOC:TP or TOC:TN (hypothesis 3), TOC-chl-a curves were estimated for each season using data subdivided into quartiles of SVA 420 , TOC:TP and TOC:TN. TOC:TP and TOC:TN were expressed as molar ratios; TOC:TP and TOC:TN data were log 10 -transformed before averaging or conducting statistical tests to avoid biased results from the inherently asymmetric properties of ratio data (Isles 2020).
To estimate temporal trends in variables of interest, TOC, log 10 (TOC:TN), and log 10 (TOC:TP) were first normalized to z-scores by (1) calculating annual average concentrations for each lake, and (2) subtracting the mean and dividing by the standard deviation of those averages for each lake (c.f. Isles and others 2018a). Mean z-scores for each site for each season and year were plotted, allowing an assessment of relative changes in lake variables as well as relative differences between seasons. Loess regressions were used to fit trends in the data, with error bars representing 95% confidence intervals around the estimate of the mean calculated as 1.96 standard errors.
To predict spatially explicit changes in chl-a concentrations with a 1.7 mg L -1 increase in TOC concentration, regression models of chl-a as a function of TOC, TOC 2 (to account for the unimodal response of chl-a to TOC), SVA 420 , and log 10 (TOC:TP) or log 10 (TOC:TN) were fit using the LTM data and the linear model (lm) function in R. The best-performing and most parsimonious models were used to predict changes in chl-a.
The CLIML lakes were divided into 47 regions using k-means clustering based on their geographic coordinates (coordinate reference system SWEREF 99). The clustering produced uneven numbers of sites per region due to differences in the density of sampled lakes in different regions, ranging from 12 lakes on the island of Gotland, to 227 lakes in a densely sampled area in southern Sweden, with a median of 103 lakes per region. Median values were calculated for TOC, SVA 420 , log 10 (TOC:TN), and log 10 (TOC:TP) in each region. The chl-a given a 1.7 mg L -1 increase in TOC was then calculated based on the regression models described above (we assumed that TOC:TP remained unchanged).
Regression, multiple regression, correlation, and cluster analyses were conducted using R (R Core Team 2018).

Shape of the TOC-chl-a Curve in Different Seasons (Hypothesis 1)
A quadratic relationship was the best fit to the chl-a and TOC data when all seasons were aggregated (AIC = 38,951, 39,026, and 40,153 for quadratic, Gaussian, and linear fits, respectively); however, few of the observations had TOC greater than that at the apex of the curve (Figure 1). The Gaussian curve was the next-best fit, and the linear fit was the weakest. The maximum value of the quadratic fit occurred at TOC = 28 mg L -1 , and the maxi-mum value of the Gaussian fit occurred at TOC = 19 mg L -1 , suggesting that the choice of functional form may strongly affect the estimated position of the maxima. The quadratic and Gaussian curves were both more similar to the binned medians than the linear curve ( Figure 1). The quadratic curve was also the best fit in all individual seasons , with Gaussian curves providing the next-best fit in summer and autumn, and a linear curve providing the next-best fit in spring. Because quadratic curves were identified as the best fits for all seasons, quadratic curves were also fitted seasonal TOC-chl-a curves subsetted to quartiles of the TOC:TP ratio, TOC:TN ratio, and SVA 420 .
Although there was a significant, but relatively weak unimodal relationship observed for the LTM dataset as a whole, a unimodal relationship was more-clearly observed in the fall season , where chl-a declined quickly at TOC above 20 mg L -1 , but not in the spring and summer, when chl-a concentrations increased with TOC across the range of observed data. This may partially reflect that volumetric rather than areal chl-a concentrations were used. Using our estimates of areal chl-a concentrations for the summer season, the TOC-chl-a relationship changed from a monotonic increase to a unimodal relationship and Gaussian (red solid line) curves fit to the data. Models were fit to the median in non-log space, but fits are presented here on log scaled y-axis. Medians of equally spaced bins along the TOC gradient are also shown (black circles) for bins containing at least 10 observations. Vertical lines at each median represent the 25th and 75th percentiles of chl-a for each bin. Figure 2. TOC plotted against chl-a in the LTM lakes and sampled lakes partitioned by season, absorbance, and stoichiometry. We hypothesized that there would be a unimodal relationship, with the height of the curve controlled by carbon-nutrient stoichiometry, and the horizontal position of the maximum controlled by the light climate (seasonality and TOC color). Panels A-C: Seasonality of the TOC-chl-a relationship in the LTM lakes. Individual observations are small gray circles, with median values of equally spaced bins along the TOC axis shown as the large circles and error bars representing the interquartile range. Blue dashed lines indicate best-fitting curves to data from each season (in all cases quadratic models). Blue diamonds indicate the maxima of fitted curves. Panels D-F: Summer TOC plotted against chl-a from the LTM lakes, with colors of points representing quartiles of TOC:TN (panel d) and TOC:TP (panel e), and fall TOC plotted against chl-a showing quartiles of SVA 420 . (panel f). Dark blue represents the lowest quartile of each variable, light blue is the second quartile, light red the third, and dark red is the highest quartile of each variable. Large circles represent medians of equally spaced bins along the TOC axis for each subset of the data; bins with less than 5 samples were discarded. Dashed lines represent quadratic regressions fit to each subset of the data. with a maximum at 19.9 mg L -1 TOC (Figure 3). However, most observations (88%) remained below the maximum of the curve.

Light Absorption Effects on TOC-chl-a Relationships (Hypothesis 2)
Unimodal relationships were only observed in the highest bin of SVA 420 , and therefore, we could not identify the effect of SVA 420 on the horizontal position of the maximum chl-a. However, SVA 420 did not appear to be related to the shape of the TOC-chl-a curve ( Figure 2F, Supplemental Figure 2), given that the fitted curves and binned medians were quite similar in all quartiles. Lower quartiles (< 50%) of SVA 420 were only marginally higher than upper quartiles (> 50%) of SVA 420 at low TOC concentrations. In the upper quartile of SVA, the chl-a peak occurred at a TOC concentration of 15 mg L -1 .

Stoichiometry Effects on TOC-chl-a Relationships (Hypothesis 3)
Stoichiometry was related to the height of the TOCchl-a curve ( Figure 2D,E Supplemental Figure 2). In summer, when chl-a was highest and nutrient limitation was likely to be strongest, lower quartiles of log 10 (TOC:TN) and log 10 (TOC:TP) were associated with higher chl-a concentrations at a given concentration of TOC ( Figure 2D, E). Log 10 (-TOC:TP) appeared to have a stronger effect on the shape of the TOC-chl-a curve than log 10 (TOC:TN), which was similar but less pronounced.
TOC:TP stoichiometry also had a strong effect on the areal estimates of summer biomass (Figure 3). When areal biomass was binned by quartiles of TOC:TP stoichiometry, only the lowest quartile of TOC:TP showed a monotonic increase, whereas the upper three quartiles reached maxima within the range of the observed data. The unimodal character of the relationship was most important in the highest quartile of TOC:TP, where a large fraction (28%) of the observations had TOC concentrations above the maximum. In all other quartiles, fewer than 1% of observations had TOC concentrations above the modeled peak of the unimodal curve.

Spatial Patterns in TOC, Absorption, and Stoichiometry
TOC, SVA 420 , log 10 (TOC:TN), and log 10 (TOC:TP) were lowest in the northern subarctic and mountainous regions (Figure 4). Log 10 (TOC:TN) and log 10 (TOC:TP) were extremely low in the southernmost tip of Sweden where the landscape is predominantly agricultural. Log 10 (TOC:TN) was also relatively low in most of southern Sweden, where atmospheric N deposition is highest. SVA 420 has generally similar patterns to log 10 (TOC:TP) and log 10 (TOC:TN), with low values in the northwest subarctic areas and generally high values in central and southern Sweden, with the exception of the southernmost region and the area around Stockholm in the southeast (Figure 4).
TOC concentration was positively correlated with SVA 420 (r = 0.57, p < 0.0001) in the CLIML dataset, with SVA 420 increasing with TOC up to a TOC concentration of approximately 20 mg L -1 , then remaining approximately constant with further increases in TOC. However, there was no consistent correlation of TOC concentration with SVA 420 over time within individual lakes from the LTM dataset (Supplemental Figure 4), suggesting a disconnect between spatial patterns and temporal changes.
Temporal Trends in TOC, Absorption, Stoichiometry and chl-a TOC increased over the monitoring period, with periodic cycles at an approximately decadal scale around a general increasing trend ( Figure 5). Log 10 (TOC:TN) decreased between 2007 and 2018 and was consistently lower in spring relative to summer and autumn. Log 10 (TOC:TP) increased with an apparent step-shift toward higher ratios around 2006-2008. Chl-a had a generally increasing trend from 1996 to 2014 during summer, but this trend declined from 2014 to 2018. There were no trends apparent in chl-a during the spring or the fall. There was not a clear link between temporal changes in summer chl-a and temporal changes in TOC, TOC:TN or TOC:TP across all regions (Figure 5). However, the mean within-lake correlation of TOC and chl-a was significantly positive during summer (Supplemental Figure 5). There was not a clear spatial pattern in sites with increasing chl-a, with the largest increases and decreases both observed in southern Sweden (where sites were more likely to be impacted by anthropogenic nutrient stress) (Supplemental Figure 6).

Predicted Changes in chl-a with Continued Browning
The best-performing models based on the LTM data for all seasons predicted log 10 (chl-a) as a function of TOC, TOC 2 , and log 10 (TOC:TP). The models explained more variance in summer than in spring or fall (Table 1). Adding SVA 420 did not improve model predictions for any season, and the fitted coefficients for SVA 420 were not significant. Using log 10 (TOC:TN) instead of log 10 (TOC:TP) reduced the R 2 of the models by 0.09, 0.08, and 0.11 for spring, summer, and fall, respectively. The best performing models for each season (Table 1, Supplemental Figure 7) were applied to the regional median values of TOC and log 10 (TOC:TP) (Figure 4) to predict log chl-a concentration in each region. The models were then used to predict chl-a with an increase of 1.7 mg L -1 in TOC (the mean change in TOC in the LTM data across all lakes for the period of this study) and then to calculate the percent change in chl-a for each region and for each season (Figure 6). Percent changes ranged from -4.48 to + 37.8% in the spring (in absolute units Dchl-a ranged from -0.23 to + 1.5 lg/L chla, median 0.52), from + 4.1 to + 76.6% in summer (0.30-2.5 lg/L chl-a, median 0.82), and from -1.12 to + 30.0% ( -0.06 to + 1.34 lg/L chl-a, median 0.46) in fall. Proportional changes in chl-a were highly variable spatially ( Figure 6), with the largest increases observed in the northeastern regions, moderate increases in central Sweden, and only modest increases or decreases in most of southern Sweden where TOC was highest (Figures 4 and 6).

DISCUSSION
There was a unimodal relationship between TOC and chl-a concentrations under some conditions, but under most conditions chl-a increased with b Figure 5. Temporal changes in TOC, SVA 420 , log 10 (TOC:TN), log 10 (TOC:TP), and chl-a from 1996 to 2018. All variables were converted to z-scores within each lake to compare relative changes across lakes. Loess smoothing lines are for data subsetted to each season (spring = April-June, summer = July-August, fall = September-October). Shaded areas around each loess fit represent 95% confidence intervals around the mean value. All model coefficients were highly significant (p < 0.0001). Figure 6. Percent changes in chl-a (colors) with a 1.7 mg L -1 increase in TOC for each region calculated based on median regional values of TOC and log 10 (TOC:TP) from the CLIML dataset (see Figure 3).
increasing TOC (Figure 1). The unimodal response was less consistent than we expected based on previous studies (Karlsson and others 2015; Kelly and others 2018; Bergströ m and Karlsson 2019), but our finding is supported by a recent analysis that found increases in PP with increasing DOC across large gradients in boreal Canada (Bogard and others 2020), by a recent study which identified no unimodal relationship between DOC and chl-a in a survey of lakes in central Ontario (Senar and others 2021), and by targeted experimental studies which commonly observed increases in GPP with increasing DOC (Klug 2002;Zwart and others 2016). A unimodal fit was significantly better than a linear fit between TOC and log 10 (chl-a) for the dataset as a whole (Fig 1), and for all seasons individually (Figure 2), but there was substantial variability around the fitted lines. A unimodal relationship between TOC and chl-a concentrations was mostly observed in nutrient poor lakes with high SVA 420 and high carbon-tonutrient ratios (Figures 2, 3). There was a strong positive response of chl-a to increasing TOC at low concentrations of TOC (below 5-10 mg L -1 ) in all seasons. In the spring and summer, there was a positive relationship between TOC and chl-a even at high TOC concentrations, suggesting that the cotransport of nutrients with DOC resulted in a fertilization effect under most circumstances (Klug 2002;Kissman and others 2016). The strongest unimodal response was observed in the fall (when light availability to pelagic phytoplankton are low due to low surface irradiance and deeper mixed layers), and this response was restricted to the 50% of lakes with the highest SVA 420 , TOC:TN and TOC:TP. This partially supports hypothesis 1 (the TOC concentration at which the peak occurs varies seasonally); however, rather than affecting the location of the peak, seasonal differences affected whether a peak occurred at all.
The unimodal relationship was somewhat stronger when considering areal (depth-integrated) rather than volumetric chl-a concentrations (Figure 3). Increasing DOC reduces the mixed layer depth (Fee and others 1996;Kelly and others 2018), and therefore even if volumetric chl-a concentrations show modest increases with DOC, declines in areal chl-a may be observed, resulting in a stronger unimodal shape. Indeed, this effect is central to the current theories explaining the unimodal relationship between DOC and PP (Klug 2002;Kelly and others 2018), and when we applied our estimate of areal chl-a we observed a unimodal curve for some quartiles of TOC:TP that was not present in the volumetric estimates. This distinction between volumetric and areal concentrations may help to reconcile our results with several other studies, which have revealed a stronger unimodal response of GPP to increasing DOC (Hanson and others 2008; Ask and others 2012; Seekell and others 2015b). However, even after controlling for estimated effects of DOC on thermocline depth, there were very few samples with DOC higher than the peak of the fitted unimodal curve, suggesting that fertilization is the dominant effect across most conditions. This is consistent with previous work. We observed the peak of the unimodal curve for summer areal chl-a at 19.9 mg TOCAEL -1 ; for context, Sobek and others (2007) estimated that only about 4.6% of lakes globally had DOC higher than 20 mgAEL -1 . It is possible that our models were somewhat biased due to the use of TOC as a proxy for DOC. In particular, some samples might have contained high phytoplankton biomass (and thus TOC and chl-a) or other particulate organic matter even though DOC was relatively low, which could result in some points in the upper right quadrants of the panels on Figures 1, 2 and 3, potentially skewing the fitted curves and shifting the unimodal peaks to higher levels of TOC. This may partially explain the increase in binned medians of chl-a at high TOC concentrations in some seasons (for example, Figure 2, panel e). However, in general the binned medians aligned very well with the fitted curves (Figure 2), suggesting that this effect was probably modest.
The height and the horizontal position of the peak of the TOC-chl-a curve were not clearly explained by SVA 420 , because maxima were not observed for most quartiles of SVA 420 and because there was no clear separation among the fitted curves or the binned medians for the four quartiles ( Figure 2, Supplemental Figure 2). We therefore conclude that hypothesis 2 (the TOC concentration at which the peak chl-a occurs varies with SVA 420 , with higher SVA 420 resulting in maximum chl-a at lower DOC concentrations), was not supported. Although our analysis did not identify a relationship between SVA 420 and the shape of the chl-a response curve, the importance of SVA 420 may be masked by underlying regional differences in landcover (Supplemental Figure 1) and spatial correlations with DOC stoichiometry (Figure 4) which have stronger influence on phytoplankton biomass. The effect may also have been masked by confounding effects of lake size and lake residence time. Lake size is expected to impact light availability (and thus the horizontal position of the maximum of the DOC-chl-a curve; compare Kelly and others 2018), because larger lakes have deeper Lake Phytoplankton Respond to DOC Concentration and Nutrient Stoichiometry mixed layers than smaller lakes (due to longer fetch and greater response to wind mixing), and consequently lower average irradiance in the mixed layer. Larger lakes with longer residence time may also be expected to have more highly processed (and therefore less colored) DOC than smaller lakes (Kothawala and others 2014). Experimental studies explicitly testing the role of DOC with different absorbance characteristics would help to resolve this question.
The height of the TOC-chl-a curve was explained in large part by carbon-nutrient stoichiometry, supporting hypothesis 3(that is, the height of the curves was controlled by stoichiometry, but a peak was not always evident). The clear vertical separation among curves binned by both TOC:TN and TOC:TP demonstrated that the fertilization response to dissolved organic matter loading is affected by carbon-to-nutrient stoichiometry. For both nutrients, lower carbon-to-nutrient ratios were associated with a steeper proportional response of chl-a to increasing TOC, and higher carbon-to-nutrient ratios were more often associated with unimodal responses (Figure 2, panel e, Supplemental Figure 2). These responses were similar to those reported in earlier field studies (Jones 1992;Klug 2002;Deininger and others 2017), and support previous studies, suggesting that light limitation, and corresponding unimodal responses to browning, may be particularly common in nutrient-poor northern lakes (Karlsson and others 2009, Supplemental Figure 3). Even after estimating areal chl-a, a unimodal relationship was mostly observed in the highest quartiles of TOC:TP during summer ( Figure 3). This result suggests that relatively higher nutrient levels may allow positive net growth rates of phytoplankton even in low-light conditions, whereas relatively nutrient-poor conditions lead to suppression of growth rates by both light and nutrient limitation.
The stronger effect of TOC:TP compared to TOC:TN on the response of chl-a to browning likely reflects the fact that TN is a poor indicator of bioavailable N in Swedish lakes due to the high proportion of biologically unavailable organic N ( Morris and Lewis 1988;Bergströ m 2010;Berggren and others 2014), whereas a larger proportion of the TP pool is bioavailable. A better proxy for truly bioavailable N might therefore exhibit a stronger relationship with trends in chl-a, as dissolved inorganic nitrogen earlier has been shown dominate over phosphate as the limiting nutrient for phytoplankton in nutrient poor boreal lakes (Elser and others 2009;Bergströ m and others 2015;Isles andothers 2018, 2020).
There were strong spatial patterns in all of the factors proposed to influence the shape of the TOCchl-a curve (SVA 420 , TOC:TP, and TOC:TN), and these patterns showed consistent regional differences in Sweden (Figure 4), reflecting the diversity of landscape characteristics across the large climatic, geological, and vegetational gradients in the study area. The observed patterns have implications for predicting the strength of the phytoplankton response in different regions to further increases in DOC. The mountainous and subarctic regions in the northwest of Sweden had low values for all of the measured variables. Low log 10 (-TOC:TP) and log 10 (TOC:TN) ratios in the north were likely due to differences in landcover between the sub-arctic biome (which tends to be dominated by shrub birch forests and tundra heath vegetation) and the boreal forests (spruce and pine) which dominate most of central Sweden (Supplemental Figure 1) and are associated with high DOC loads (Sobek and others 2007). Log 10 (TOC:TP) and log 10 (TOC:TN) were both relatively high in central and eastern parts of northern Sweden, reflecting the high prevalence of both boreal forest and wetlands in these regions. In the south, there were contrasting patterns of log 10 (TOC:TP) and log 10 (-TOC:TN). Log 10 (TOC:TN) was relatively low in most of southern Sweden, which is probably due to high N deposition in this region (compare Isles and others 2018) and higher agricultural inputs of nutrients (Sponseller and others 2014). However, log 10 (TOC:TP) was relatively high in most of the south (excluding the area surrounding Stockholm in the southeast likely impacted by urbanization, and Skå ne in the southernmost tip of Sweden which is heavily agricultural).
We observed substantial long-term variation in log 10 (TOC:TN)and log 10 (TOC:TP) ratios which, based on our analysis, we expect to have substantial impact on lake productivity. These temporal patterns in TOC-nutrient stoichiometry, which have not to our knowledge been previously shown, would have an important role in mediating responses of lake productivity to browning. The trends in log 10 (TOC:TN) and log 10 (TOC:TP) did not match trends in TOC, indicating that other factors were differentially impacting TN and TP concentrations. SVA 420 also did not match trends in TOC, indicating that, with respect to both stoichiometry and color, the factors that influence DOC properties are not the same as those that influence DOC quantity, as has recently been observed elsewhere (Jane and others 2017). Given the strong correlation of carbon-nutrient stoichiometry with the chla response to TOC (Figure 2, Table 1, Supplemental Figure 2), it was expected that long-term changes in log 10 (TOC:TP) and log 10 (TOC:TN) would have influenced lake phytoplankton biomass. There was no obvious relationship between temporal changes in carbon-nutrient stoichiometry and chl-a. However, chl-a is influenced by a number of factors not considered here (for example, temperature, atmospheric or direct anthropogenic nutrient deposition uncorrelated with DOC loads, and food web structure), which may obscure any underlying influence of browning and nutrient stoichiometry.
Future browning (if it occurs) would likely lead to increases in epilimnion chl-a concentrations in most lakes, with the greatest increases observed during summer. Our models predicted that an increase of 1.7 mg L -1 TOC (the median change over the study period) could increase chl-a by up to 76% in some areas in northern Sweden, and by a substantial fraction (30-60%) in much of central Sweden ( Figure 6). Somewhat smaller increases would likely be expected if we considered areal rather than epilimnetic chl-a concentrations (Figure 3); unfortunately, lake area data were not available for the CLIML lakes and it was additionally not clear how many lakes in the dataset stratified annually, so we were not able to do this analysis. Projected increases in chl-a concentrations were somewhat lower in the fall than in the summer in these regions, possibly reflecting greater light limitation during these months (Figure 2c). By contrast, lakes in southern Sweden tended to have high TOC, high SVA 420 , and high TOC:TP, and therefore would likely experience only modest increases or decreases in epilimnetic chl-a in response to future browning.
It is likely that other northern regions experiencing browning (Monteith and others 2007) would also generally experience increases in epilimnetic chl-a with continued browning, particularly because many other areas have lower median DOC concentrations than Sweden (Sobek and others 2007), so are more likely to be on the increasing side of the DOC-chl-a curve. It is also true that some northern regions are not experiencing browning, particularly in coastal areas where high and increasing precipitation may result in dilution of lake DOC (de Wit and others 2016a, b). In Sweden, previous studies have indicated the largest increases in DOC in the south, whereas some northwestern areas have had minimal changes or even declines (Isles and others 2018). Our findings are most relevant to lakes whose catchments have minimal human development (urbanization or agriculture), because the effects of cotransport of nutrients with DOC are likely less important when there are other, direct sources of nutrients added to lakes. Our results have the most significant implications in boreal and subarctic areas, where a large fraction of the world's lakes are located (Verpoorter and others 2014), direct human impacts on catchments are often minimal, and climate change is likely to result in wide-ranging changes to DOC concentrations and composition (Creed and others 2018).
Our results clearly indicate that unimodal relationships between DOC and chl-a are possible and likely generally present, but that most lakes are on the left side of this curve, and will experience increasing chl-a concentrations with continued browning. It is unclear whether these patterns reflect changes in PP across DOC gradients. Several studies of pelagic PP have observed maximum production at lower DOC concentrations than the maxima identified in our study for chl-a, even after accounting for areal chl-a concentrations (Hanson and others 2003;Seekell and others 2015a). PP and phytoplankton biomass are related but not synonymous, and it has been observed that the highest rates of PP may be observed at intermediate chl-a concentrations (Carpenter and Kitchell 1984). However, across large gradients, chl-a and PP have been found to be well correlated, with R 2 as high as 0.8 (Smith 1979), and a recent study examining PP across over 170 lakes in boreal Quebé c found monotonic increases in PP across DOC gradients from around 1-30 mg C L -1 (Bogard and others 2020), consistent with our findings. Bergströ m and Karlsson (2019) also observed peaks in epilimnetic phytoplankton biomass at lower DOC concentrations than observed in this study; however, most lakes in this study were quite brown and nutrient poor, so this result is consistent with what we observed in the higher quartiles of TOC:TP (Figures 2,  3), which reached peak chl-a at lower TOC than lakes with lower TOC:TP. Experimental manipulations testing theories of the trade-off of light and nutrient dynamics are complicated due to the difficulty and expense of inducing changes in lake DOC concentrations at enough levels to resolve a unimodal peak, particularly when other covariates like TOC:TP are examined; however, future studies should attempt to better clarify the relationship Lake Phytoplankton Respond to DOC Concentration and Nutrient Stoichiometry between chl-a and PP across DOC gradients, to determine the degree to which our findings can be used to make conclusions regarding lake primary productivity.
Although our findings suggest that most lakes will have increased epilimnetic chl-a if browning continues, further research into the mechanisms driving this response is needed. In particular, our model ( Figure 6) assumes constant DOC-nutrient stoichiometry over time, an assumption which is unlikely to be met given the trends in TOC:TN and TOC:TP ( Figure 5). However, the apparent contrast in trends in TOC:TN and TOC:TP suggests that changing DOC-nutrient stoichiometry could either increase or decrease median chl-a concentrations, making it unclear whether our estimates should be revised upwards or downwards. Additionally, cotransport of nutrients with dissolved organic matter, which is at the core of our analysis, is just one mechanism by which browning affects lake processes. Our model should be combined with models integrating the effects of thermal structure to better understand changes in areal chl-a concentrations across a range of lake mixing regimes, as well as models of the effects of temperature on lake phytoplankton productivity and temperature-driven changes in phytoplankton stoichiometric requirements, particularly because browning may be changing the identity of limiting nutrients in lakes toward increased N limitation (Isles and others 2020). Future studies attempting to quantify phytoplankton responses to browning should also consider the effects of browning on food web structure and top-down controls (Kelly and others 2016; Creed and others 2018). mits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit h ttp://creativecommons.org/licenses/by/4.0/.

D A T A AV AI L ABI LI TY
All public monitoring data are available for download from the Swedish long-term monitoring program (http://miljodata.slu.se/mvm/).