Physical characteristics of northern forested lakes predict sensitivity to climate change

Evidence suggests that boreal-lake ecosystems are changing rapidly, but with variable ecological responses, due to climate warming. Paleolimnological analysis of 27 undeveloped northern forested lakes showed significant and potentially climate-mediated shifts in diatom communities and increased carbon and biogenic silica burial. We hypothesize the sensitivity of northern forested lakes to climate change will vary along two physical gradients: one reflecting direct, in-lake climate effects (propensity to thermally stratify), the other reflecting indirect watershed effects (watershed to lake-surface area ratio). We focus on the historical response of algal communities to test our two-dimensional sensitivity framework. Historical algal response was summarized by measures of diatom community turnover, changes in species and diagnostic species groups, and measures of siliceous algal and overall primary production (biogenic silica, carbon burial). Measures of algal production increased across all lake types, with carbon burial proportionately higher in polymictic lakes. Greater diatom community change occurred in deep, stratified lakes with smaller watersheds, whereas diatom species groups showed variable responses along our two-dimensional sensitivity framework. Physical characteristics of lakes and watersheds could serve as predictors of sensitivity to climate change based on paleo-indicators that are mechanistically linked to direct and indirect limnological effects of climate change.


Introduction
The boreal biome is the world's second largest ecoregion. It covers approximately 15% of the earth's land surface and contains > 60% of the world's fresh surface waters in over 3 million lakes (Ruckstuhl et al., 2008;Schindler & Lee, 2010;Verpoorter et al., 2014). Based on multiple lines of evidence, northern forested lake ecosystems are changing rapidly, with unprecedented appearances of cyanobacterial blooms (Winter et al., 2011;Cottingham et al., 2015;Favot et al., 2019), significant shifts in algal communities (Rühland et al., 2008(Rühland et al., , 2015Edlund et al., 2011;Saros et al., 2012;Hobbs et al., 2016), increased watershed inputs of dissolved solids (organic carbon and nitrogen: Stottlemyer & Toczydlowski, 2006;Keller, 2007;Zhang et al., 2010;Corman et al., 2018), changes in thermal structure (Richardson et al., 2017), and increased carbon burial in lake sediments (Anderson et al., 2013;Heathcote et al., 2015;Bartosiewicz et al., 2019). Given the lack of local human impact on most boreal lakes, climate change and atmospheric deposition of nitrogen, phosphorus and other pollutants are the most likely drivers of these ecological shifts. However, while regional-scale drivers such as climate may be expected to generate moderately coherent change in ecosystems, prior work on northern forested lake ecosystems has highlighted the high degree of spatial and temporal variability in the ecological response (i.e., sensitivity) of lakes to climate change over the past century. For example, paleolimnological data spanning the twentieth century reveal that some lakes have experienced more than a 50% turnover in diatom species assemblage starting around 1920 (Saros et al., 2012), but with no apparent change in total production, while other lakes from the same region have experienced an increase in algal production and a shift to cyanobacterial blooms starting in the 1970s Edlund et al., 2017). This variability over different scales is common in many regions, and is often qualitatively attributed to watershed and/or within lake differences (Baron & Caine, 2000;Benson et al., 2000;Velle et al., 2005;Patoine & Leavitt, 2006;Richardson et al., 2017).
Many studies have ascribed ecological changes in boreal lakes to climate change through correlative means (e.g., coherence of trends among lake-sediment records; Rühland et al., 2008), but with less consideration of the mechanisms involved. Climate may affect lakes through direct effects on ice cover (Jensen et al., 2007), surface temperature (Schneider & Hook, 2010;Dokulil, 2014;Richardson et al., 2017), thermal structure (Kraemer et al., 2015;Richardson et al., 2017), wind-induced turbulence (Austin & Colman, 2007), and mixing Richardson et al., 2017), which in turn can alter light regime and nutrient regeneration from stratified bottom waters (Keller, 2007). Autotrophic communities and primary production are highly responsive to changes in light and nutrients in both the short-term (seasonal succession; Winder et al., 2009;Berger et al., 2010) and the long-term (decadal assemblage shifts; Rühland et al., 2008;Saros et al., 2012). Possible indirect (watershed) effects include increased inputs of nutrients (Stoddard et al., 2016) and dissolved organic carbon (DOC; Corman et al., 2018) as a consequence of increased soil microbial activity and changes in subsurface drainage in response to decreased winter snowpack and warmer temperatures (Stottlemyer & Toczydlowski, 2006). Indeed, widespread increases in watershed DOC export reported from North America and Europe have been attributed to a warming climate as well as other anthropogenic drivers (Clark et al., 2010;Larsen et al., 2011;Solomon et al., 2015;Williamson et al., 2015;Meyer-Jacob et al., 2019). Because lake response to climate drivers is clearly multidimensional, we are only beginning to be able to predict which lakes will respond most strongly to climate change and through what mechanisms (e.g., Edlund et al., 2017;Richardson et al., 2017;Bartosiewicz et al., 2019).
To determine the potential sensitivity of lakes to climate change, with sensitivity defined by degree of change in climate-linked paleoecological indicators, we propose a two-gradient ordination ( Fig. 1) based on (1) the relative influence of watershed effects, represented by the ratio of watershed to lake-surface area (WALA; stronger watershed influence on lakes with WALA > 10), and (2) lake thermal structure, represented by lake geometric ratio (GR = surface area 0.25 /maximum depth), a metric for the propensity of a lake to thermally stratify (GR > 5 polymictic, GR < 3 strongly stratified; Gorham & Boyce, 1989). We suggest that the key proximate factors (e.g., internal and external nutrient loading) linking climate change to lake response will be largely colinear with these two axes, thus providing a simplified yet mechanistic framework for understanding ecological change in northern forested lakes. The GR captures lake thermal structure, including length of ice-free season, and depth of stratification, i.e., direct climate effects; thermal structure further influences light regime and in-lake nutrient cycling (van de Waal et al., 2009). Similarly, WALA is a well-established predictor of nutrient and DOC loading as well as hydrologic residence time-key watershed-mediated physical and chemical variables in the biological function of lakes that reflect indirect climate effects. The variables of this framework, including lake depth and surface area, have also been recently highlighted as the best global predictors of the magnitude of lake stratification changes resulting from climate change (Kraemer et al., 2015;Richardson et al., 2017).
To explore the effects of climate change on boreal lakes and the utility of this framework we examine paleolimnological records from 27 lakes located in the boreal and northern forested region of the northcentral US. All lakes are in protected watersheds (tribal or national and state parks and forests) with no or minimal development in their watersheds and minimal change in land use during the last 50 years. Lead-210 dated sediment cores from each lake were examined for measurements of ecosystem function including autotrophic production, carbon burial, diatom community turnover, and species response over the last 120-150 years. Ecosystem change metrics were projected onto our direct and indirect climate response ordination to determine if sensitivity of lakes to direct and indirect climate drivers is predictable and mechanistically sound.  (Fig. 2). Lakes were characterized by standard physical and biogeochemical/limnological parameters including measures of water quality (total P, chlorophyll-a, Secchi depth, total N, nitrate-nitrite, ammonium, DOC, pH, conductivity, alkalinity, dissolved SiO 2 ) and physical measures (mean and maximum depth, geometric ratio, WALA) mostly as part of a long-term monitoring effort . Temperature records from 1890 to 2015 were compiled from meteorological stations in northern Minnesota (Eveleth, Minnesota), northern Upper Peninsula Michigan (Marquette, Michigan), and northern Lower Peninsula Michigan (Hart, Michigan). Long-term temperature records were all downloaded from the National Oceanic and Atmospheric Administration (NOAA) Global Historical Climatology Network, version 3 (http:// www. ncdc. noaa. gov/ ghcnm/ v3. php), which archives monitored records from weather stations across the United States. Data were summarized as winter (DJF; December-February mean), summer (JJA; June-August mean) and mean annual temperature (metANN).

Sample collection and analysis
Single sediment cores (0.5-2.0 m long) were collected from the central area of the main basin of each lake. Cores were collected using a drive-rod piston corer equipped with a 6.5 cm diameter polycarbonate barrel (Wright, 1991) except Siskiwit and Desor lakes, which were sampled using an HTH gravity corer (Renberg & Hansson, 2008). All cores were sectioned in the field in 0.5 to 1-cm increments; core samples were returned to the laboratory and stored at 4°C for further analyses. Freeze-dried samples of all core sections are archived at the St. Croix Watershed Research Station (Minnesota, USA).
Lead-210 was measured in all cores by isotopedilution alpha spectrometry, and dates and sedimentation rates were determined according to the CRS (constant rate of supply) model (Appleby and Oldfield, 1978). Dating and sedimentation errors were determined by first-order propagation of counting uncertainty (Binford, 1990). In each core, between 15 and 19 core sections were analyzed for Pb 210 activity to determine age and sediment accumulation rate for the past 150 years. Most dating models have been previously published (Table S1) or are shown in Fig.  S6.
Biogenic silica (BSi), the siliceous algae (diatoms and chrysophytes) content of sediments, was measured on freeze-dried subsamples (~ 30 mg) from each core section that were digested using 40 ml of 1% (w/v) Na 2 CO 3 solution heated at 85°C in a reciprocating water bath for five hours (DeMaster, 1979;Conley & Schelske, 2001). A 0.5-g aliquot of supernatant was removed from each sample at 3, 4, and 5 h. After cooling and neutralization with 4.5 g of 0.021 N HCl solution, dissolved silica was measured colorimetrically on a Unity Scientific SmartChem 170 discrete analyzer as molybdate reactive silica (SmartChem, 2012). Measured BSi concentrations were converted to flux using dry-mass sedimentation rates in each core.
Diatom samples were prepared by placing approximately 0.25 cm 3 of homogenized sediment in a 50-cm 3 polycarbonate centrifuge tube and adding 2-5 drops of 10% v/v HCl solution to dissolve carbonates. Organic material was subsequently oxidized by adding 10 ml of 30% hydrogen peroxide and heating for 3 h in an 85°C water bath. After cooling, the samples were rinsed with distilled deionized water to remove oxidation byproducts. Aliquots of the remaining material, containing the diatoms, were dried onto 22 × 22 mm #1 coverglasses, which were then permanently attached to microscope slides using Zrax mounting medium (MicrAP, Pittsburgh, USA). Diatoms were identified along random transects to the lowest taxonomic level under 1250X magnification with oil immersion optics. A minimum of 400 valves was counted in each sample. Species abundances represent percent abundance relative to total diatom counts. Identification of diatoms used regional   floras (e.g., Patrick & Reimer, 1966, 1975Reavie & Smol, 1998;Camburn & Charles, 2000) and primary literature.
Carbon burial was calculated as the product of the percent organic carbon and the Pb 210 derived drymass accumulation rate (DMAR) from each sediment core corrected for sediment focusing. Percent organic carbon was estimated by percent loss-on-ignition at 550°C multiplied by a constant (0.469; Dean, 1974;Heathcote et al., 2015). Focus corrections were applied to each core based on the ratio of the Pb 210 flux measured in the sediment to the regional atmospheric Pb 210 flux (Engstrom & Rose, 2013;Heathcote and Downing, 2012;Heathcote et al., 2015). This correction provides a more representative estimate of the basin-wide carbon burial. Finally, core sections less than 10 years in age (based on Pb 210 age) were excluded from the analysis to remove artefacts due to incomplete carbon mineralization (Gälman et al., 2008;Heathcote et al., 2015).

Calculating metrics of change
Four metrics were calculated to determine historical changes in productivity and biology in each lake. Productivity was summarized as organic carbon burial or flux (Heathcote et al., 2015), a measure that combines in lake productivity and watershed sources of carbon to lakes, and biogenic silica flux, a measure of historical diatom and chrysophyte productivity. Biological change in lakes was summarized using squared chord distance, a standardized measure of diatom community change (Bennion et al., 2004;Simpson et al., 2005). Biological change within well-known diatom indicator species and species groups was also calculated. All metrics of change were considered in relation to two physical measures of lake susceptibility to climate change-direct climate effects on lake thermal structure summarized as lake geometric ratio (GR; Gorham and Boyce, 1989) and indirect climate effects mediated by watershed processes summarized as WALA (Wetzel, 2001). Analyses were performed in R (R Core Team, 2016).
Changes in organic C flux and BSi flux at core sites were averaged over two time periods representing pre-Euro-American settlement and modern conditions (i.e., pre-1900 and post-1980; landuse histories summarized in Edlund et al., 2011). For a given parameter, the ratio of the average value of all post-1980 was calculated against the average measured values dated prior to 1900 from each core. Therefore, a value < 1 represented a decrease in the parameter over the time period, and a value > 1 an increase. Ratios of change in organic C and BSi flux were summarized on a biplot of WALA (on a log 10 scale) versus GR and scaled using bubble plots.
Floristic change in each lake's diatom community was determined between similar time periods as above using the squared chord distance (Bennion et al., 2004;Simpson et al., 2005). Squared chord distance is a dissimilarity measure with values ranging from 0 (samples are exactly the same) to 2 (samples are completely different). The squared chord distance was calculated between two time periods in each core: pre-1900 (1880 ± 12 years) to the core top or present. Chord distances were summarized on a biplot of WALA (on a log 10 scale) versus GR and scaled using bubble plots.
To explore the response of indicator diatom species and species groups, the difference between pre-1900 percent abundance and post-1980 percent abundance was calculated for the following diatom species (or groups of species): Asterionella formosa Hass., all Aulacoseira species, Fragilaria crotonensis Kitton, large cyclotelloid species, and small cyclotelloid species. The small cyclotelloids category groups small centric diatoms (typically with diameters < ~ 15 µm) once placed in the genus Cyclotella, but now considered in the genera Discostella and Lindavia/Pantocsekiella including D. stelligera   (Table 1)  These diatom species and species groups have been identified as indicators of climate-and anthropogeniclinked physical, chemical, and biological change in lakes (e.g., Winder et al., 2009;Rühland et al., 2008;Saros et al., 2012). Changes were summarized on a biplot of WALA (on a log 10 scale) versus GR and scaled using bubble plots with blue bubbles showing increases and red bubbles showing decreases in species abundances. All diatom data have been previously published or are available in the Neotoma Paleoecology Database (Table S1).

Lakes and climate
The 27 lakes encompass a broad range of water quality and physical limnological characteristics (Table 1). Lakes were generally oligotrophic to mesotrophic with total phosphorus (TP) ranging from 3 to 23 µg l −1 , chlorophyll from 0.7 to 5.5 µg l −1 , and Secchi depths of 0.75 to 8.5 m. Increasing TP was correlated with lakes having higher WALA and higher geometric ratios (Fig. S1). Lakes also varied in color from clear to strongly stained as measured by DOC, which ranged from 3.9 to 20.4 mg l −1 , with higher DOC lakes correlated with higher geometric ratios (Fig. S1). Lakes were generally circumneutral (pH from 6.4 to 8.6), weakly to well-buffered (alkalinity 5 to 160 mg l −1 ), with specific conductivity from very low to moderate levels (5 to 160 µS cm −1 ). Lake and watershed size varied across three orders of magnitude with lake areas from 10-12,000 ha and watershed areas of 70-88,000 ha. As such, WALA also varied from very low (WALA 1.5, Mukooda Lake, VOYA) to large systems such as Lac La Croix in the BWCAW with WALA of 72 (Fig. 1). Finally, the geometric ratio of the lakes varied from those unlikely to stratify with high GR (15.5, Outer Lagoon, APIS) to deep lakes with steep-sided basins that were strongly stratified, e.g., Cruiser Lake, VOYA has GR 0.9 (Fig. 1).
Regional climate records varied across our study area. Mean annual temperature increased over the period of record  in northern Minnesota and upper Michigan, but increased only after the 1970s in northern lower Michigan. (Fig. 3). Mean winter temperatures increased after the 1970s in northern Minnesota and lower Michigan, but showed little long-term trend in upper Michigan; however, all three stations show more recent temperature extremes and interannual variability in winter temperatures since the 1970s (Fig. 3). Summer mean temperatures increased in northern Minnesota from the 1890s to the 1930s and have remained relatively constant since, with notable higher interannual variability since the 1980s. In contrast, summer temperatures in upper Michigan increased steadily between 1920 and 1970, with much faster warming since the 1970s. In northern lower Michigan, a slight cooling trend between the 1930s and 1960s quickly shifted to rapid summer warming since the 1970s.

Carbon burial
Focus-corrected organic carbon burial followed two patterns of change among the study lakes (Fig. 4A). First, almost all lakes showed increased carbon burial over the last 150 years with only the deep lakes Mukooda, Trout, and Lake La Croix, and shallow Beaver Lake declining in organic C burial. Second, the lakes with the highest overall rates of C burial (2to threefold higher than the great majority of lakes) are mostly shallower lakes that rarely stratify or have weaker stratification with GR > 3.

Siliceous algal productivity
Historical diatom and chrysophyte productivity was summarized as changes in flux of BSi (SiO 2 flux; Fig. 4B) to the sediments. BSi flux increased in all study lakes with no strong pattern in magnitude of change within our two-gradient ordination. Lakes with the greatest increase in BSi flux included Cruiser, Manitou, Bass, Ek, and Swamp lakes, which are distributed in two states and three different park/ tribal management units.

Diatom community change
Diatom community change measured as chord distance from pre-1900 to present varied broadly among lakes (Fig. 5A). There was a noticeable shallow lakes also showed large community shifts (e.g., Shell, Beaver, and Bass lakes).

Diatom species and species groups
Percent abundance change in five diatom species and species groups from pre-1900 to present showed varying patterns of increase and decrease among species and study lakes. Aulacoseira species generally decreased in abundance across our study lakes (Fig. 5B). Large declines in Aulacoseira abundance were noted in Mukooda, Peary and Shoepack lakes from VOYA and in Desor and Richie lakes from ISRO. Notable increases in Aulacoseira abundance were seen in Manitou, Brown, Ryan, and Beaver lakes, which are in three different national park units.
In comparison to our sensitivity ordination, lakes with the largest increases in Aulacoseira abundance are arrayed between WALA of 4 to 10 and GR 2 to 5, whereas lakes with decreasing Aulacoseira abundance are generally more strongly stratified (GR < 4).
The abundance of small cyclotelloid species increased in all but six of the study lakes between pre-1900 and post-1980s (Fig. 5C). Lakes that showed dramatic increases in small cyclotelloid species, namely Desor (ISRO), Beaver (PIRO) and Manitou and Bass (SLBE), vary from weakly to strongly stratified, and are arrayed across the range of WALA. It is notable that small increases in cyclotelloid abundance were also commonly seen in the shallow lake systems that rarely stratify. Lakes that decreased in small cyclotelloid abundance were clustered together at very low GR < 2 and at WALA < 3. Although the relationships were not significant, changes in small cyclotelloid abundance showed some correlation with other modern water quality and physical lake characteristics including lake depth (negative) and TP (positive) (Fig. S2).
Large cyclotelloid species also increased in abundance among study lakes, with only six lakes having decreases in abundance (Fig. 5D). Siskiwit Lake, a large oligotrophic lake on ISRO, dominates the pattern of change with dramatically increased abundance between pre-1900 and post-1980s (Fig. 5D). Other lakes with increased abundance tend to be clustered among the deep, strongly stratified systems (GR < 3) and those with low WALA < 10). The few lakes that have decreased abundance of large cyclotelloids have somewhat higher WALA (WALA > 6). It should be noted that many of the lakes either have little change in large cyclotelloid species or this species group is absent from the lakes, especially in high GR and shallow systems. The change in large cyclotelloid species abundance was also negatively correlated to DOC and positively correlated to lake depth (Z max ; Fig. S3).
Where present, the abundance of Asterionella formosa increased in all but three of our study lakes (Fig. 5E). A. formosa abundance decreased in Elk, Locator, and Brown lakes, which are all located in northern Minnesota. This group of lakes is reasonably deep with strong stratification (GR < 3) and has WALA of approximately 5.0-6.5. Lakes with the greatest increase in A. formosa abundance included VOYA's Cruiser and Mukooda lakes. Other lakes with increasing A. formosa abundance are arrayed throughout the WALA scale. It should be noted that many shallow lakes and weakly stratified lakes did not have A. formosa (e.g., Wallace, Harvey, Shell, Swamp, Ryan, Speckled Trout, Outer, and Desor lakes).
The diatom Fragilaria crotonensis has a limited distribution and was found in only 10 out of our 27 study lakes. It is not found in lakes with < 10 µg l −1 TP, nor in lakes with Z mean > 12 m, and is not present in most polymictic lakes (GR > 4; Fig. S4). Where found, the abundance of F. crotonensis generally increased in the study lakes including large increases in Grand Sable (PIRO) and Manitou (SLBE) lakes, Richie and Ahmik lakes at ISRO, and Lac La Croix in Minnesota's BWCAW, which are all lakes with WALA > 3 and lakes that are generally more productive (Fig. S4). The few lakes that lost abundance of F. crotonensis were mostly stratified systems with a low WALA (e.g., VOYA's Mukooda and ISRO's Desor lakes).

Discussion
Long-term climate response of lakes is mediated by lake morphometry and hydrologic setting (Fritz, 2008). To predict the sensitivity of lakes to climate drivers we distilled potential effects of climate on lakes into a simple bivariate matrix. Direct climate effect is based on energy inputs to lakes (Leavitt et al., 2009) and is predicted by GR (Gorham & Boyce, 1989), the propensity to stratify that is calculated using known predictors of lake-climate response: lake-surface area and depth (Kraemer et al., 2015;Richardson et al., 2017). Indirect climate effects involve mass transfer to lakes that are mediated by the watershed, and is predicted by WALA (Fraterrigo and Downing, 2008). Indirect climate impacts, such as changes in lake clarity, watershed hydrology and connectivity , and nutrient and DOC loading have been linked to summer precipitation (McCullough et al., 2019) and winter severity (Stottlemyer & Toczydlowski, 2006), which are indirect lake-climate connections captured by WALA. Limnologists have benefitted from similar simplified predictive models because of their ability to cut through data and site variability, find broad-scale patterns, build predictive power, and drive actionable management (e.g., Carlson, 1977).
Paleolimnological records have long-been used to understand or reconstruct lake-climate response based on proxies that represent watershed, physical, and lake processes (Cohen, 2003;Fritz, 2008). While there are direct impacts of climate on lakes (warming, ice out, mixing), most paleolimnological proxies represent the endpoint of a much more complex response to climate and non-climate factors. When working on individual lakes, factors that are not directly controlled by climate such as watershed processes, food webs, landscape position, and lake morphometry (Webster et al., 1996;Kraemer et al., 2015;Wigdahl-Perry et al., 2016) must be considered when inferring climate impacts on lakes. A strategy to better understand climate impacts on lakes is to synthesize large numbers of records to highlight regional or global patterns of response to climate drivers, whether with monitored data (Schindler et al., 1996;Kraemer et al., 2015;Richardson et al., 2017) or paleolimnological records (Laird et al., 2003;Rühland et al., 2008;Heathcote et al., 2015) as we have done using paleoindicators of primary production, community composition, and species changes. In our case, the response of each paleoindicator across our bivariate matrix predicts and can be linked mechanistically to how northern forested lakes may respond to direct and indirect climate effects.
Within our 27 remote study lakes, paleolimnological proxies responded in largely predictable ways along our direct (GR) and indirect (WALA) climate susceptibility gradients. As such, we can predict not whether lakes will respond to climate, but how they might respond to direct and indirect climate drivers. For example, siliceous algal productivity has increased in all study lakes across our matrix, i.e., a generalized lake response to both direct and indirect climate drivers. Carbon burial similarly increased in almost all study lakes; however, there is an enhanced direct climate effect in weakly stratified or polymictic lakes (GR > 3), which showed larger increases in C burial. Diatom community change was greatest in strongly stratified (GR < 3) and smaller WALA (WALA < 10) lakes indicating a strong direct climate response and weaker indirect climate response. In contrast to productivity or community measures, indicator species or group response varied across the bivariate matrix. Small and large cyclotelloids responded strongly, but in opposite direction, to direct climate effects and weakly to indirect effects with strongly stratified and small WALA lakes showing decreased small cyclotelloids and increased large cyclotelloids. Changes in Aulacoseira and Asterionella formosa were less predictable along either axis; however, most lakes showed declines in Aulacoseira and increases in A. formosa abundance. Response of Fragilaria crotonensis was better predicted by WALA and not GR, indicating its utility for detecting indirect climate effects.
Organic carbon burial to lake sediments provides a historical record of changes in in lake productivity and watershed sources of organic carbon (Heathcote et al., 2015). In our study lakes, most sites showed increased rates of carbon burial between pre-1900 and recent times. Only four lakes showed decreased rates of carbon burial and these changes were all small. Lakes that showed the greatest changes in carbon burial (2-to fourfold increases) had GR > 3 and were shallow or weakly stratified. A similar negative correlation between increased carbon burial and lake depth was noted across Minnesota lakes (Dietz et al., 2015). Importantly, in undeveloped watersheds, rates of carbon burial also reflect significant negative latitude and positive temperature relations within the Anthropocene (Heathcote et al., 2015). Significant increases in epilimnetic temperatures after 1960 (especially in shallow lakes; Edlund et al., 2017) coupled with browning of lakes, appear to be interacting to promote widespread increased burial of organic carbon in northern wilderness lakes (Bartosiewicz et al., 2019).
In many lakes, anthropogenic land-use disturbance is the primary driver of 4-to fivefold increases in organic carbon burial as watershed nutrient inputs enhance lake eutrophication through agricultural intensification and erosion and transport of terrestrial carbon (Heathcote & Downing, 2012;Anderson et al., 2013;Dietz et al., 2015;Heathcote et al., 2014). Links to climate change, latitude-temperature, or basin morphometry are weak in comparison to signals of anthropogenic disturbance. However, in northern regions that have not undergone extensive landscape modification, organic carbon accumulation rates still show marked twentieth century increases-that may be linked to historical logging and land clearance (Dietz et al., 2015; but see Laird & Cumming, 2001)-but are also associated with climate warming and atmospheric deposition of reactive nitrogen (Heathcote et al., 2015).
Siliceous algal productivity, summarized as flux of BSi, increased in all of our study lakes between pre-1900 and post-1980. In anthropogenically impacted regions, increased BSi flux is readily linked to historical point source nutrient loading (Conley et al., 1993;Edlund et al., 2009) and land-use changes that enhance nutrient delivery (e.g., agriculture; Heathcote et al., 2014). Fewer studies have considered BSi flux in remote regions and its potential relationship to climate impacts. Previous reports that include some of our study lakes (Edlund et al., 2017) echo the overall patterns of increasing BSi flux we show in remote lakes. Three mechanisms likely drive these increases in productivity. Modeled and measured records suggest that length of ice cover is decreasing resulting in longer algal growing seasons (Jensen et al., 2007), epilimnetic temperatures are increasing in many remote and worldwide lakes (Schneider & Hook, 2010;Hadley et al., 2014;Edlund et al., 2017;Bartosiewicz et al., 2019) leading to enhanced algal productivity (Schaum et al., 2017), and lake and watershed responses to warming climate may shift nutrient loading and in-lake dynamics to promote productivity (Stottlemyer & Toczydlowski, 2006;Horn et al., 2011;Maier et al., 2018;Corman et al., 2018). It should be noted that warming temperatures should not dramatically increase BSi mineralization or alter the record of BSi flux as is possible with carbon burial (Bartosiewicz et al., 2019).
Large-scale and rapid shifts in algal communities are commonplace in response to changing land-use and new nutrient sources. Climate impacts on algal communities can be equally dramatic , but are often subtler and require long-term monitoring or multi-proxy paleoecological analysis to account for other anthropogenic stressors (Rühland et al., 2015). Among our study lakes, the deep, stratified lakes with small watersheds showed the greatest diatom community change since the late 1800s. This community shift encompasses not only recent periods of climate change, but also, depending on individual lakes, some level of landscape change that has variously included logging, land clearance, mining, cottage development, and reforestation that were a consequence of Euro-American settlement (Edlund et al., 2011). To understand changes since the 1800s, Hobbs et al. (2016) identified timing of significant shifts in community turnover among national parks (including 20 of our 27 lakes). At ISRO and SLBE, significant community turnover occurred in the early 1900s likely in response to early landscape changes noted above. In addition, all park units showed a significant shift as some point from the 1950s (SLBE, APIS, PIRO) through the 1970s-1990s (VOYA, ISRO) . An initial analysis of post-1970s diatom community changes in our lakes confirms broad change in communities since the 1970s in all lake types, but especially those with smaller WALA < 10 (compare Figs. 5A and S5). With all of the study lakes under land-use protection by the 1970s, these more recent changes capture periods of significant increases in mean annual air temperature (Fig. 3) and direct climate limnological impacts such as modeled increased epilimnetic temperatures and interannual variability in some of our lakes (VOYA and ISRO; Edlund et al., 2017). Sorvari et al. (2002) similarly showed diatom community changes were temporally consistent among regional lake groups in Finland, but that the type of community shift differed across lake and catchment types, indicating the need to consider species and species group changes to better determine the mechanisms driving change.
Diatom ecological groups, species complexes, and individual species have been used in experimental, mesocosm, and paleolimnological efforts to decipher climate impacts on lakes. Sorvari et al. (2002) showed the historical shift between benthic and planktonic diatom communities in Finnish arctic lakes was most strongly correlated to spring warming trends over the last 150 years. Shifts in planktonic composition to greater abundance of cyclotelloid species at the expense of larger Aulacoseira species have also been linked to limnological responses (e.g., longer ice-free season, water column stability) associated with warming winters in some northern lakes (Rühland et al., 2008(Rühland et al., , 2015. Saros and Anderson (2015) synthesized the resource requirements based on experimental and mesocosm studies of cyclotelloid species to understand the mechanisms behind species' modern and paleolimnological responses to climate-linked limnological changes in temperature, light, and nutrients. We took a similar approach to these studies using a synthesis of paleolimnological records to disentangle how five diatom taxa (Fragilaria crotonensis, Asterionella formosa) or species complexes (Aulacoseira species, small cyclotelloids, large cyclotelloids) responded to direct and indirect climate impacts as defined by our sensitivity framework.
The genus Aulacoseira is common in freshwaters worldwide, inhabiting both alkaline and low pH water bodies (Spaulding & Edlund, 2008). Species environmental preference varies from polymictic and eutrophic (e.g., A. granulata), to mesotrophic conditions (e.g., A. ambigua, A. subarctica), and includes low-alkalinity forms (e.g., A. distans). Their heavily silicified cell walls join these taxa together as a species group that requires a well-mixed water column, limiting their seasonality to periods of spring and fall turnover or shallow and polymictic conditions (Lund, 1954). There was an overall decrease in abundance of Aulacoseira taxa across most study lakes over the last century. Only four lakes showed notable increases in Aulacoseira abundance, and those four were arrayed between WALA 4 to 10 and GR of 2 to 5, suggesting that neither the propensity for stratification nor watershed size was a strong predictor of increased abundance. Rather the decrease in Aulacoseira abundance across many lakes has been used to link limnological change to climate warming in temperate lakes, particularly shorter duration of ice cover and weakened mixing due to stronger or more rapid stratification (Sorvari et al., 2002;Rühland et al., 2015).
Long-term monitoring provides one link between climate drivers and response of Aulacoseira taxa in temperate lakes. Horn et al. (2011) linked a 35-year record of plankton, meteorology, and physicochemistry in a German reservoir to show that the abundance of A. subarctica (O.Müll.) E.Y.Haw. was greater in years with short mild winters, shorter winter ice cover, earlier ice out, and especially with longer periods of spring circulation.
Paleolimnological connections with monitored limnological change and climate drivers show Aulacoseira abundance in cores negatively correlated with mean annual temperature and positively correlated with earlier ice out in Canadian temperate to tundra lakes Rühland et al., 2008).
Other paleolimnological records highlight the variability of Aulacoseira response among lakes types. Sorvari et al. (2002) studied five Finnish lakes experiencing similar climate change, with two of the lakes showing marked decreases in Aulacoseira abundance linked to meteorological and tree ring records of warming. The other three lakes showed little change in Aulacoseira abundance or did not have Aulacoseira taxa in their floras. Boeff et al. (2016) in a study of sediment cores and modern limnology from three lakes in Maine (USA), concluded that a decrease in Aulacoseira abundance in two of the lakes was indicative of stronger stratification based on modeled response. However, in a third lake, a seasonal separation of modern abundance patterns in Aulacoseira (typically a spring form) and Discostella (a summer and fall form) was not clear and the paleolimnological record showed little evidence of change in Aulacoseira abundance or inferred stratification change. The key in the Maine lakes was understanding the seasonality of taxa in the individual lakes; outlier lakes with regard to Aulacoseira autecology did not follow the widely held convention that decreased Aulacoseira abundance indicates earlier ice out and longer/ stronger stratification.
Although we do not have clear understanding of seasonal Aulacoseira ecology in our study lakes, the widespread decrease in Aulacoseira abundance supports regional limnological responses to warmer winters, earlier ice out, and stronger stratification. Edlund et al. (2017) reported significant winter warming in this same region between the 1960s and present and used thermal lake modeling to show that lakes in VOYA and ISRO had warmer summer epilimnetic temperatures and more frequent and longer thermoclines in recent years, especially in deeper lakes. Earlier ice out has been documented in the region as well (Jensen et al., 2007), although our study lakes typically do not have extensive ice out records. The decrease in Aulacoseira in our study lakes is among the most consistent metrics of change, but is clearly not a universal response among boreal lakes. The modern ecology of lakes that show increased abundance of Aulacoseira should be studied carefully to determine why they are outliers in this response (Boeff et al., 2016).
The small cyclotelloids comprise a group of centric diatom species earlier placed in the genus Cyclotella; however, more recent systematic treatments place these diatom taxa within the genera Discostella and Lindavia/Pantocsekiella . The taxa are often limited in distribution to oligoand mesotrophic lakes (Saros & Anderson, 2015), and are well represented in the deep chlorophyll layer community, especially in large lakes (Fahnenstiel & Glime, 1983). Their small size may limit their sinking rates and provide them with higher nutrient uptake (greater surface area to volume ratio; SA: V) and photosynthetic capacity, allowing greater growth rate (Saros and Anderson, 2015;Rühland et al., 2015;Reavie et al., 2017). Across our study lakes, small cyclotelloids increase in 21 of 27 lakes between pre-1900 and recent, including many shallow lakes with GR > 4. Those that did not show increases are mostly deep, oligotrophic, and well stratified lakes with small WALA.
The increase in small cyclotelloid species across northern hemisphere lakes has generally been linked to direct climate effects including decreased ice coverage, increased growing season length, stronger thermal stratification (Sorvari et al., 2002;Rühland et al., 2003Rühland et al., , 2008Winder et al., 2009;Enache et al., 2011), and longer turnover (Wiltse et al., 2016). However, among our study lakes, increases in small cyclotelloid abundance occur along both our direct climate gradient (GR) and our indirect watershed gradient (WALA). And systems that we might expect to be strongly and directly affected by climate (deep, strongly stratified, low WALA), in fact show decreased abundance of small cyclotelloid species, suggesting that a more nuanced consideration of which species are actually responding and in what lakes may be warranted in future studies.
A similar conclusion was reached by Saros and Anderson (2015) and Reavie et al. (2017) with regard to directionality of change among lakes, speciesspecific changes, and timing of change. Overall, the timing of increases in small cyclotelloid species has been shown to have occurred earlier (late 1800s) in high latitude lakes and later (the 1970s) in temperate regions (Rühland et al., 2008(Rühland et al., , 2015; however, other systems have exhibited later or earlier shifts (Rühland et al., 2015;Reavie et al., 2017). In the Laurentian Great Lakes, changes in cyclotelloids are overlain with varying histories of other anthropogenic impacts and recovery. However, positive correlation between cyclotelloid abundance and meteorological data such as minimum annual temperature suggests a mechanistic, limnological linkage between cyclotelloid abundance and direct climate drivers . Saros and Anderson (2015) proposed a different strategy to understand how cyclotelloid abundances might be interpreted in the context of global change. They argue that the autecology of individual species needs to be considered, including seasonality, habitat preference, nutrient requirements and limitation, and light and DOC response. Saros et al. (2014), Malik and Saros (2016), and Malik et al. (2017) studied the temperature, light, and nutrient requirements of arctic populations of Lindavia radiosa, L. bodanica (both large cyclotelloids in our analysis), L. comensis, L. ocellata, and Discostella stelligera and found large differences in resource interactions among species. Lindavia ocellata was a low light taxon, L. comensis responded to nutrient addition except at high light and temperature, L. radiosa did best with nutrients under high light, L. bodanica did well with nutrient addition regardless of temperature, and Discostella stelligera did well in moderate light under P-limitation, but when under N/P co-limitation required low light. Malik et al. (2018) found that timing of ice-off could explain interannual differences in abundance of L. bodanica and D. stelligera in a Maine (USA) lake, but that the diatoms were responding predictably throughout the ice-free period to thermal structure, light condition, and nutrient limitation. Processes that alter light availability (including water clarity and water column stability) and nutrient concentrations likely play a major role in controlling the growth of small centric diatoms in lakes. Similar differences in resource ecology and seasonality have shed light on the variability of cyclotelloid response in the Great Lakes  and UK lakes (Thackeray et al. 2008). Coupling neo-limnological and experimental autecology with paleolimnological approaches provides a more robust and mechanistic interpretation of potential climate-driven limnological and biological responses (inferred mixing depth; Saros et al., 2012;Boeff et al., 2016;Strock et al. 2019) across lake types and contrasting climate conditions.
Although traditionally placed in the genus Cyclotella, the large cyclotelloid species are currently treated in the genus Lindavia and comprise taxa such as L. radiosa, L. bodanica, L. comta, L. lemanensis, and L. intermedia (Spaulding et al., 2019). These taxa are broadly distributed in stratified oligo-to mesotrophic lakes (Reavie & Kireta, 2015;Saros & Anderson, 2015), where they often form a characteristic community within the deep chlorophyll layer (Malik et al., 2018) and control their sinking rate through the production of extensive fibrils of β-chitin and external mucilage (Herth & Barthlott, 1979). The abundance of large cyclotelloid species increased slightly in most of our study lakes where this group was present (in 15 out of 21 lakes), with greatest increases in deep, strongly stratified lakes with lower WALA, especially in Siskiwit Lake (ISRO). The six lakes that showed declines were from among three national park units and had generally higher WALA and were not strongly stratified. Based on mesocosm studies, large cyclotelloid species exhibit strong interactive responses to light, temperature, and nutrients (Saros et al., 2012;Malik & Saros, 2016;Malik et al., 2017). Although there is some variability among species (Malik et al., 2017), the large cyclotelloids generally respond positively to nutrient addition and have greater growth rates in meta-and hypolimnetic waters (Saros & Anderson, 2015).
Other studies have reported similar recent increases in the abundance of large cyclotelloid species. Saros et al. (2012) and Strock et al. (2019) modeled changes in thermocline depth based on dramatic increases in large cyclotelloid species in Siskiwit and Desor lakes at ISRO, inferring that it had deepened 6 m and 3 m, respectively, in response to greater wind speed in the twentieth century (Austin & Colman, 2007). Similar increases of large cyclotelloids in many of our study lakes may be related to wind speed in or near the Great Lakes. However, increases also occurred at inland lakes in northern Minnesota (Little Trout, Mukooda, and Trout) where wind speed has moderated during the twentieth century (Edlund et al., 2017). As our study is the first to consider these changes in cyclotelloid abundance among many lake types and regions, it is important to recognize that other factors, notably a stable thermocline and perhaps increased external and internal nutrient loading may also contribute to change, and that wind speed alone cannot be solely responsible for the widespread increase in large cyclotelloids in the northern lakes. In contrast, most Great Lakes show decreased abundance of large cyclotelloids linked to warming surface temperatures and intensified stratification, although multiple anthropogenic stressors confound the diatom response .
The colonial araphid diatom, Fragilaria crotonensis, is widely distributed in the plankton of meso-or eutrophic lakes (Morales et al., 2013) where it forms band colonies of 10 to 50 cells imparting some resistance to sinking (Rühland et al., 2015). In our data set and elsewhere (Bradbury, 1988), F. crotonensis is absent from small, shallow lakes, as well as from deep lakes (Z mean > 12 m) with low nutrients. Across our bivariate sensitivity matrix, Fragilaria crotonensis occurred in only 10 of 27 lakes. Where found, it increased in abundance in lakes with WALA > 3, and decreased in abundance mostly in lakes with smaller WALA, suggesting that indirect watershed effects are driving changes in this species. The connection of climate to indirect watershed effects is likely multifaceted. Across the northern hemisphere, an increase in atmospheric deposition of reactive nitrogen (Nr) (Galloway et al., 2008;Holtgrieve et al., 2011;Baron et al., 2013) has driven changes in remote ecosystems such as boreal and alpine lakes (Saros et al., 2005;Heathcote et al., 2015;Spaulding et al., 2015). Although most monitoring stations from our study area show a decrease in Nr deposition since 2000 , the response to elevated Nr inputs is likely longer-term and chronic in nature. Stottlemyer and Toczydlowski (2006) determined the indirect mechanism linking climate change and watershed export of Nr in these northern ecosystems based on long-term watershed-scale monitoring. They argue that elevated Nr deposition has increased the watershed burden of Nr above historical (natural) background and that warming winters with lower precipitation has increased mineralization and export of DIN, DON, and DOC to receiving waters. Modern studies on the ecology of F. crotonensis indicate that this diatom does well during spring mixing when P: Si nutrient ratios are low (Bradbury, 1988) and is thus receptive to N additions especially under low P availability (Saros et al., 2005;Michel et al., 2006). In our northern forested lakes, especially mesotrophic systems, winter warming provides the link between climate and watershed export of N and C (Stottlemyer and Toczydlowski, 2006).
Asterionella formosa is a widespread araphid species found in the plankton of oligo-to eutrophic lakes throughout the northern hemisphere (Lund, 1949;Michel et al., 2006). Its long, narrow cells grow in characteristic star-shaped colonies that impart a high SA: V for enhanced nutrient uptake and a low sinking rate (Rühland et al., 2015). In many lakes it is a common spring bloom-former (Lund 1949), but can persist in the plankton well into the summer months (Malik et al., 2018), where it is vulnerable to fungal parasitism and grazing (Winder and Schindler, 2004;Kagami et al., 2007).
Among all diatom taxa or species groups within our dataset, Asterionella formosa showed the most consistent pattern of increased abundance in lakes (where it was present). Asterionella formosa was absent from many shallow and weakly stratified lakes. Only three lakes showed declines in A. formosa, and these were deep, stratified, N-limited mesotrophic systems (Table 1; Bergström, 2010) with a WALA of 5-6. The widespread increase in abundance across our sensitivity matrix suggests that A. formosa is responding to both direct and indirect environmental drivers. In other regions increases in A. formosa appear to be a direct result of atmospheric Nr deposition (Saros et al., 2005;Wolfe et al., 2006;Spaulding et al., 2015) or eutrophication (Horn et al., 2011;Thackeray et al., 2008). However, increases in A. formosa abundance have also been noted under conditions of declining nutrient levels (Sivarajah et al., 2016), raising questions as to the possible role of climate in controlling the growth response of Asterionella.
Reported climate responses by Asterionella formosa are confounding. Malik et al. (2018) studied a lake in Maine (USA) over two years that differed by late and early ice out, but stratified on nearly the same date. Asterionella abundance was much higher in the late ice out year, but cells persisted in the water column longer in the early ice-off year, perhaps indicating the importance of light or limiting nutrients (silica, P). In the UK Lake District, Thackeray et al. (2008) showed that the spring maximum of A. formosa occurred earlier with warming. Sivarajah et al. (2016) suggested that water column warming, longer periods of open water, and changes in mixing might have contributed to the increase in A. formosa in an Ontario lake in the absence of nutrient addition.
Indirect climate drivers are also identified in controlling abundance of A. formosa in Swedish lakes. Maier et al. (2018) used a 12-year sediment trap sequence and showed that A. formosa increased in the sediment record in years characterized by high runoff events, i.e., indirect climate effects, that preceded fall turnover.
Among lakes supporting A. formosa populations, P enrichment from point sources, watershed sources, or rapid destratification can lead to abundant A. formosa (Thackeray et al., 2008;Maier et al., 2018;Horn et al., 2011). In oligotrophic lakes, especially N-limited systems, A. formosa responds equally well to N addition from anthropogenic sources (Saros et al., 2005) and climate-linked sources such as glacial meltwaters (Slemmons et al., 2015(Slemmons et al., , 2017 and watershed export (Stottlemyer and Toczydlowski, 2006). The response of Asterionella formosa in our study lakes likely reflects a complex response to nutrient availability controlled by both direct and indirect climate effects. Our dataset includes both P-and N-limited systems, such that resolving the response of this species will require lake-specific limnological studies (e.g., Thackeray et al., 2008;Saros et al., 2012;Malik et al., 2018 Lakes are environmental sentinels that can help us detect, monitor, and understand the threat of climate change to our aquatic resources (Williamson et al., 2009). Most efforts to understand climate sensitivity of lakes have been guided by comparing meteorological and limnological monitoring records (Kraemer et al., 2015;Richardson et al., 2017;McCullough et al., 2019); these records are often limited in temporal extent, geographic extent, or discrepancies in monitored parameters among sites. Paleolimnological records extend the geographic and temporal range of studies and offer common analytical strategies to measure multiple limnological responses and assess the sensitivity of lakes to direct and indirect climate effects. Paleolimnological approaches generate multiple lines of evidence that provide ample opportunity to define lake susceptibility to direct and indirect climate drivers. We show that diatom productivity measured as BSi flux is increasing in all lakes across both direct and indirect climate effect gradients, whereas increases in carbon burial among remote lakes are more prevalent in polymictic systems. Diatom community changes most strongly reflect direct climate impacts as deep and stratified remote lakes show greatest community turnover. Diatom species or species groups vary in their response across our direct/ indirect climate matrix. For example, Aulacoseira species and large cyclotelloids responded strongly along the direct climate gradient with decreasing abundance of Aulacoseira and increased abundance of large cyclotelloids in the strongly stratified lakes. Other species such as Fragilaria crotonensis increased in abundance in lake with high WALA, a response to indirect climate drivers, whereas Asterionella formosa increased in abundance across our susceptibility matrix, a positive response to both direct and indirect climate drivers. We welcome broader testing of our susceptibility matrix and response variables. With northern forested regions holding over 60% of the world's fresh surface water, understanding the sensitivity of these lakes to climate change is paramount. funded by the US National Park Service (Cooperative Agreement P17AC01133), the Great Lakes Network Inventory and Monitoring of the US National Park Service under Cooperative Agreement P20AC00364-01, the Air Resources Division of the US National Park Service (Agreement: H6000082000), the US Geological Survey (PM15-94151), the Environment and Natural Resources Trust Fund through the Minnesota Department of Natural Resources (Subcontract B33299), and the Grand Portage Band of Lake Superior Chippewa through the US Geological Survey (8607-AZJ00). Sadly, co-author Dr. Norman Andresen passed away in August 2019 while this manuscript was being written.
Funding This study was funded by the Great Lakes Network Inventory and Monitoring of the US National Park Service under Cooperative Agreement P20AC00364-01, the Air Resources Division of the US National Park Service (Agreement: H6000082000), the US Geological Survey (PM15-94151), the Environment and Natural Resources Trust Fund through the Minnesota Department of Natural Resources (Subcontract B33299), and the Grand Portage Band of Lake Superior Chippewa through the US Geological Survey (8607-AZJ00).

Data availability
The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.

Declarations
Conflict of interest Co-Author Jasmine E. Saros is an Associate Editor of Hydrobiologia.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits 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 http:// creat iveco mmons. org/ licen ses/ by/4. 0/.