Inconsistent browning of northeastern U.S. lakes despite increased precipitation and recovery from acidi ﬁ cation

. Multiple studies have reported widespread browning of Northern Hemisphere lakes. Most examples are from boreal lakes that have experienced limited human in ﬂ uence, and browning has alterna-tively been attributed to changes in atmospheric deposition, climate, and land use. To determine the extent and possible causes of browning across a more geographically diverse region, we examined watercolor and dissolved organic carbon (DOC) time series in hundreds of northeastern U.S. lakes. The majority of lakes have increased in both DOC and color, but there were neither coherent spatial patterns in trends nor relationships with previously reported drivers. Color trends were more variable than DOC trends, and DOC and color trends were not strongly correlated, indicating a cause other than or in addition to increased loading of terrestrial carbon. Browning may be pronounced in regions where climate and atmospheric deposition are dominant drivers but muted in more human-dominated landscapes with a limited extent of organic soils where other disturbances predominate.


INTRODUCTION
Many freshwater ecosystems in the Northern Hemisphere have been browning, that is, becoming enriched in chromophoric dissolved organic carbon (DOC) in response to a combination of changing atmospheric deposition, climate, and land use (Monteith et al. 2007, Clark et al. 2010, Kritzberg 2017. Darker watercolor can have strong effects on several lake attributes including carbon (C) budgets, light penetration, primary productivity, stratification, and thermal regimes (reviewed by Solomon et al. 2015). Browning trends have been widely documented in the Northern Hemisphere, with the strength and causes of these trends varying geographically (Monteith et al. 2007, Clark et al. 2010, Finstad et al. 2016. Causes of browning are inferred based on a relatively small sample of lakes representing a narrow range of geographic conditions dominated by cold climates, low population density, and organic soils. In these systems, recovery from acidification has been suggested as the primary driver in regions affected by atmospheric deposition of sulfate (Driscoll et al. 2003, Monteith et al. 2007, Haaland et al. 2010), whereas changing climate or land cover may be more important in regions that have been less affected by acid deposition (Eimers et al. 2008, Freeman et al. 2014, Kritzberg 2017. Browning trends are thus common in relatively remote northern regions, but may not be as widespread as often assumed (Filella and Rodríguez-Murillo 2014), especially in regions covering broader and more heterogeneous environmental gradients.
Browning has been reported both in terms of increasing DOC concentrations, or increasing color (measured as visible light absorbance by chromophoric DOC), often interchangeably. Recent studies have shown, however, that concentrations of DOC and color may be unrelated, especially in systems with high retention times or high productivity where colored DOC is preferentially lost and uncolored DOC produced in situ is replenished (Weyhenmeyer et al. 2012, Massicotte et al. 2017, Evans et al. 2017. In these cases, temporal trends in DOC and color may not be the same (Jane et al. 2017). Therefore, in regions where the concentrations and proportions of terrestrial DOC in lakes may be lower than in regions where browning has typically been reported, browning may be less pronounced and increasing DOC may not result in increasing color. We hypothesized that this happens in regions where human activities lead to increased nutrient loading, and hence, where lake DOC dynamics may be driven more by production and processing in the water than in the catchment.
Here, we investigate broadscale DOC and color trends in a region of the northeastern and midwestern United States where changes in climate and atmospheric deposition have occurred, but land use and other local factors also have a strong influence on lake carbon dynamics . First, we quantified the magnitude and frequency of increasing DOC and color in hundreds of lakes with a wide range of land cover and human activity, where previously identified drivers of lake browning have changed over the last 2-3 decades. Second, we compared DOC and color trends with trends in nitrogen, phosphorus, and chlorophyll a to determine whether browning and eutrophication occur simultaneously. Finally, we evaluated which lake, landscape, and climate properties were correlated with trends in DOC and watercolor.

Data source and data filtering
We used the LAke GeOSpatial Database for the Northeastern United States (LAGOS-NE v1.087.1; Soranno et al. 2015Soranno et al. , 2017 accessed through the LAGOS-NE R package (Stachelek et al. 2017) as our source of limnological and geographic data. All data details, including units and protocols for inclusion in LAGOS, are specified in Soranno et al. (2017). Monthly precipitation data to calculate trends were accessed from a supplemental climate dataset that accompanies the LAGOS-NE database . Limnological response variables are summarized in Table 1, and geospatial data used as predictors are summarized in Table 2. The data in LAGOS-NE originate from various sources and have been controlled for quality; the details of the harmonization and control procedures can be found in Soranno et al. (2015). From this database, we used annual means of DOC, watercolor (hereafter, color), total phosphorus (TP), total nitrogen (TN), nitrate-N (NO 3 ), and chlorophyll a (hereafter, chlorophyll). We used true color (measured on filtered samples) when available, but for a subset of data, only apparent color (measured on unfiltered samples) was available. Out of 56,434 color values we used, 42,283 were from true color, and 14,151 were converted from apparent color based on the linear relationship between true and apparent color when both were measured (R 2 = 0.77, slope = 0.72, intercept = −5.15). In few instances, where apparent color was reported at, or close to the detection limit (1006 out of 15,157 potential values) this estimation yielded negative true color estimates; we removed these values from our analyses. Data from the State agency category were not v www.esajournals.org included in this analysis because of low resolution in apparent color that yielded poor correlation with true color, and because true color was measured in most lakes.
To evaluate trends, we restricted the dataset to lakes that had at observations of a given limnological variable from five separate years in a time period that spanned at least 10 yr between 1980 and 2013. Lakes that met the filtering criteria for DOC and color only occurred in the northern states of the LAGOS-NE study extent (Fig. 1), which can differ strongly in environmental characteristics such as land use  and connectivity (Fergus et al. 2017) from the southern states. Hence, we further restricted the datasets for nutrients and chlorophyll to lakes that occurred in the same 8-digit hydrologic unit code (HUC 8) as lakes with data for DOC and color, resulting in comparable geographic distributions of lakes for all variables (see Appendix S1: Figs S1-S4 for maps of lakes with nutrient and chlorophyll data to compare with DOC and color data in Fig. 1).
A comparable data filtering approach has been used with the LAGOS database to evaluate temporal trends in nutrients . For that study and ours, the robustness of the analyses lies in the very large number of individual ecosystems, which compensates for the difficulty of capturing directional trends for the lakes with a more limited number of observations and temporal data coverage. Our assumption is that a strong enough, systematic influence of wellknown drivers of browning that have changed systematically over time should induce comparable patterns in browning over the population of lakes despite some variability in individual time series.
Our filtering approach resulted in data from 254 to 1714 lakes for each variable (Table 1). The data used to estimate trends spanned at least 10 yr and had observations from a median of 12 yr for all variables except TN and chlorophyll, which had slightly lower median number of years with data, with 9 and 11 yr, respectively. Very few lakes had sufficient data to calculate a Note: The median time series reflects the number of years of data included for each lake, which at a minimum needed to include observations from two years that spanned a 10-yr period.  v www.esajournals.org trend for both DOC or color and any of the nutrients, so analyses for each variable were conducted on a different population of lakes.
For predictor variables (summarized in Table 2), we used individual watershed data when they were available in LAGOS-NE-GEO. Soil carbon data were not available in LAGOS, so we used soil C data from the global soil organic carbon map from the Food and Agriculture Organization from the United Nations (FAO) soils portal (FAO 2012, Kochy et al. 2015. In LAGOS, some climate and hydrology metrics are not calculated for individual lakes, so we used HUC-12 watershed-scale data. Most HUC-12 watersheds included only one study lake, and we only used HUC-12 scale data for variables with broad spatial structure (Lapierre et al. 2018) that should not be biased by not making a calculation for individual lakes. Percent change in precipitation was calculated for the time period corresponding to the respective DOC or color trends in each lake. Sulfate deposition data are only calculated for LAGOS on a 5-year interval, so we used trends calculated for the entire time frame of the study, rather than for a different time period for each lake that corresponded to limnological data. The lack of overlap between lakes with time series data for DOC, color, and nutrients made it difficult to use nutrient trends as a predictor for DOC or color trends, but we did include median total phosphorus (a single average for each lake) as a predictor that should be related to trophic status.

Statistical methods
We used Bayesian hierarchical linear models to estimate temporal trends in lake DOC, color, nutrients, and chlorophyll (Gelman and Hill 2007). Bayesian hierarchical models were also used to estimate trends in climate variables (annual precipitation and temperature) that were then used as potential predictors of browning. These models also allowed us to evaluate potential drivers of observed spatial variability of estimated trends for DOC and color. A separate model was fitted for each variable. The model was as follows: where y i is log e -transformed DOC (or other response variable) for observation i, α j and β j are the intercepts and slopes that describe the temporal trend for lake j given the predictor of year (x i ), respectively, γ α 0 is the grand mean intercept across all lakes, γ β 0 and η p are the intercept and slopes describing the relationship between geographical drivers (Z p ) and the estimated temporal trends (β j ), respectively. Diffuse normal priors were used for γ α 0 , and γ β 0 , and a Bayesian LASSO (least absolute shrinkage and selection operator; Park andCasella 2008, Tibshirani 2011) prior was used for the slope parameters, η p , where slope parameters received a double exponential (DE) LASSO prior (η p )~DE(λ), where λ~EXP(10). All predictors Z p were standardized prior to the analysis. A diffuse uniform prior was used for the residual standard deviation (α), and we modeled the variance-covariance matrix (Σ) using the scaled inverse-Wishart distribution (Gelman and Hill 2007). We ran three parallel Markov chains beginning each chain with random starting values. Each chain was run for 10,000 iterations, from which the first 5,000 were discarded (Brooks and Gelman 1998).
We determined predictor significance by evaluating whether the 95% credible interval of the coefficient overlapped with zero. Further, to better summarize uncertainty in parameter estimates, we calculated the posterior probability that the estimated trends (β j ), and effects of predictors (η p ) on those trends, were in the same direction as the estimated posterior mean. The proportion of explained variance was calculated for each level of the following methods outlined in Gelman and Pardoe (2006). Models were fitted by calling the program JAGS (Plummer 2003) using the jagsUI package (Kellner 2018) in the program R (R v. 4.0.3 Core Team 2018).

Widespread but inconsistent browning of northeastern U.S. lakes
A majority of lakes had positive trends in both DOC and color (73% and 62%, respectively), and v www.esajournals.org 36% of lakes had a high (>90%) posterior probability of increasing in DOC and 27% of lakes had a high posterior probability of increasing in color. Besides a cluster of positive and highly probable color trends in the northeast, there were no obvious spatial patterns in the strength and direction of the trends (Fig. 1).
Not all increases in DOC are associated with true browning; that is, water color did not increase as consistently as DOC, and DOC and color were not increasing concurrently. For the 306 lakes with estimated trends for both DOC and color, there was a weak positive relationship between the trends, with only 5% of variance in DOC trends explained by color trends (Appendix S1: Fig. S7). This does not mean that DOC and color are not correlated among lakes, or over time within lakes, but rather, that most changes in DOC did not result in a corresponding change in color for the studied lakes. Just over half of lakes (59%) had DOC and color trends in the same direction, and these lakes were mostly found in a relatively narrow geographic extent around the Adirondack region toward the northeast of the study area (Fig. 1c).

More lakes are increasing in carbon than nutrients
Among the limnological variables in our analysis, DOC showed the most positive trends among the studied lakes (Fig. 2), with lakes in the northeast being more likely to have increased in DOC compared with nutrients or chlorophyll (Fig. 2). The observations that (1) DOC increased in a higher percentage of lakes than color and (2) there was no strong correlation between DOC and color trends (Appendix S1: Fig. S7) suggest that trends in DOC are not entirely driven by increased loadings of terrestrial, colored carbon. Additionally, there was no significant relationship between DOC or color trends and nutrients or chlorophyll (linear models, P > 0.05 for all comparisons), the latter explaining no more than 3% of variance in DOC or color trends (Appendix S1: Fig. S8).

Lake and landscape characteristics do not explain trends in DOC and color
We were unable to explain spatial variability in DOC and color trends using well-established drivers of browning. We related trends to several lake and landscape characteristics, including trends in climate (annual mean precipitation and air temperature) and sulfate (SO 4 ) deposition (percent change from 1980 to 2010) that have been used to explain browning trends in other studies: SO 4 deposition decreased in 100% of study lakes, annual mean precipitation increased in 97% of study lakes, and annual mean air temperature increased in 97% of study lakes. Despite these consistent shifts in drivers over the past three decades, there was no meaningful relationship between any driver and trends in either DOC (Fig. 3) or color (Fig. 4).
We evaluated the potential effects of static lake and landscape characteristics but found no significant predictors of DOC trends. More productive lakes in wetter, cooler regions with higher percent wetlands in their catchments appeared to become more colored, but only a very small amount of variance in color trends (6%) was explained by changes in annual mean precipitation, annual mean temperature, percent wetlands in the catchment, and median TP (Fig. 4). Overall, despite systematic increases in precipitation, temperature, and recovery from acidification across the study extent, those changes and other drivers were generally unrelated to changes in DOC or color. Fig. 3. Changes in DOC as a function of each of the 16 predictor variables. Solid points are the posterior mean lake-specific trend (AE95% credible interval). The solid line represents the posterior mean of the hierarchical regression line, and the shaded area is the 95% credible region. Numbers below the predictor variable name are the probability that the effect of the covariate is in the direction (positive or negative) of the posterior mean. None of the variables had posterior probability ≥0.90. AMP stands for annual mean precipitation and AMT for annual mean air temperature.

DISCUSSION
Browning of northeastern U.S. lakes is decoupled from recovery from acidification and changes in climate We found that some level of browning has occurred in hundreds of northeastern U.S. lakes over the past several decades, but unlike previous studies, the increases in concentrations of dissolved organic carbon (DOC) and watercolor were virtually unrelated to trends in climate and atmospheric deposition. This is intriguing because annual mean temperature and precipitation have systematically increased while Fig. 4. Changes in color as a function of each of the 16 predictor variables. Solid points are the posterior mean lake-specific trend (AE95% credible interval). The solid line represents the posterior mean of the hierarchical regression line, and the shaded area is the 95% credible region. Numbers below the predictor variable name are the probability that the effect of the covariate is in the direction (positive or negative) of the posterior mean. Predictor variable names and probabilities in blue represent those with a posterior probability ≥0.9. AMP stands for annual mean precipitation and AMT for annual mean air temperature. atmospheric deposition of sulfates has systematically decreased over the study area, and such trends have been associated with increasing production and mobilization of terrestrial organic carbon in soils and from land to water, respectively (Freeman et al. 2001, Monteith et al. 2007, Erlandsson et al. 2008. There is no reason to believe that these processes are not expressed in our study area, but it appears that other factors override color, and mostly DOC, dynamics in northeastern U.S. lakes. In particular, if the main source of additional DOC is not terrestrial, or if there have been changes in DOC processing in lakes over time due to changes in composition (Jane et al. 2017), DOC trends would be decoupled not only from trends in climate and atmospheric deposition, but also from trends in color.
In this study, we show such a decoupling in an area where anthropogenic sources of nutrients related to agriculture and atmospheric deposition are higher than in most previous studies that have documented browning of aquatic ecosystems. This decoupling, along with the absence of spatial patterns in DOC and color trends, suggests that watershed-scale factors linked to inlake production and processing of DOC predominate over regional-to continental-scale drivers.
What then explains the increasing concentrations of DOC and color that have been observed in the majority of studied lakes? It is possible that nutrient-driven increases in primary production led to the production of an autochthonous, noncolored source of DOC (Finlay et al. 2009, Evans et al. 2017). However, we found no relationship between nutrient or chlorophyll trends and browning trends in the small subset of lakes with concurrent time series (Appendix S1: Fig. S8), and no clear eutrophication trends in the study area that could explain an increase in internal production; overall, lake trophic changes are much less prominent in the region compared with DOC and color ( Fig. 2; Oliver et al. 2017). These results suggest that eutrophication and subsequent in situ DOC production alone cannot explain observed increases in DOC.

Reconciling browning in boreal vs. northeastern U.S. lakes
We propose two non-mutually exclusive hypotheses that could explain the different response of northeastern US vs boreal lakes to changes in climate and atmospheric deposition. First, browning has mostly been documented for lakes surrounded by organic soils (see Monteith et al. 2007, Evans et al. 2012, Kochy et al. 2015, and changes in soil chemistry may simply not result in the increased mobilization of meaningful amounts of organic carbon in carbon-poor soils. As opposed to several studies conducted in boreal Canada and Europe, there is often close to no organic soil coverage in our study extent, except at the northern edge where it rarely exceeds 10% coverage (Kochy et al. 2015). There was no significant relationship between DOC and color trends vs. soil C (Figs. 3, 4), presumably due to the small range in soil C compared with boreal landscapes (Kochy et al. 2015). However, percent wetland coverage in lake catchments was the only landscape predictors with a significant (positive) effect on color trends. Although small wetlands may not have been captured in broader estimations of soil C (Kochy et al. 2015), this result indicates that browning can be observed in our study area when large amounts of mobile organic carbon are present in surrounding landscapes.
Second, multiple interacting drivers in this region may have opposite effects on color and DOC, thereby weakening the response to the widespread and consistent changes in climate and sulfate deposition. Previously studied boreal regions are not exempt from human influence and may have experienced land-use conversions, forestry, and drainage (Evans et al. 2012, Kritzberg 2017. There is, however, little evidence that these pressures would obscure the broader scale browning trends observed in these landscapes. Several lakes in northeastern United States, on the other hand, have faced residential, agricultural, and urban developments that could simultaneously alter the loading of terrestrial DOC from land to water and the relative contribution of terrestrial vs autochthonous DOC (see Evans et al. 2017). Even in the absence of direct human effects, consistent changes in the composition of DOC have been observed in Wisconsin lakes in the northeastern United States despite little to no change in DOC concentrations (Jane et al. 2017), suggesting that the biogeochemical pathways affecting lake DOC have changed over time. Therefore, even if DOC and color in northeastern U.S. lakes responded similarly to boreal lakes to v www.esajournals.org climate and recovery from acidification, this response could be hidden among numerous and potentially interacting changes in the drivers of lake chemistry. This, in turn, could explain the inconsistent browning of northeastern U.S. lakes compared with their boreal counterparts, because the effects of climate and atmospheric deposition drivers, which operate rather homogeneously over very broad spatial scales (Lapierre et al. 2018), are muted by more local drivers related to land cover and land use.
Both hypotheses simultaneously explain our results. In particular, while trends in DOC and color in northeastern U.S. lakes may not respond as strongly to climate change and recovery from acidification compared with northern lakes, we still observed browning in the majority of the lakes. However, it is unclear where this DOC comes from, and increasing DOC rarely led to darker watercolor in the present study, and lakes that became darker (i.e., increasing color) rarely increased in DOC. Carbon cycling in lakes from our study area thus appears to differ from previous studies conducted in higher altitudes and latitudes, which have often coincided with high abundance of organic soils and minimal (or different types of) human influence beyond changing atmospheric deposition and climate. A better understanding of geographic context would help discern whether a phenomenon such as browning can be readily extrapolated across regions.

ACKNOWLEDGMENTS
We thank members of the Cross-scale interactions (CSI) Limnology and Continental Limnology groups for fruitful discussion and for developing the LAGOS-NE database. This research was supported by the Macrosystems Biology Program at the U.S. National Science Foundation (awards 1638679, 1638554, and 1638539). JFL is supported by a NSERC-Discovery grant. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. Jean-Francois Lapierre and Sarah M. Collins contributed equally to this manuscript.

LITERATURE CITED
Brooks, S. P., and A. Gelman. 1998. General methods for monitoring convergence of iterative