Landscape Controls on Nutrient Stoichiometry Regulate Lake Primary Production at the Margin of the Greenland Ice Sheet

Global change is reshaping the physical environment and altering nutrient dynamics across the Arctic. These changes can affect the structure and function of biological communities and influence important climate-related feedbacks (for example, carbon (C) sequestration) in biogeochemical processing hot spots such as lakes. To understand how these ecosystems will respond in the future, this study examined recent (< 10 y) and long-term (1000 y) shifts in autotrophic production across paraglacial environmental gradients in SW Greenland. Contemporary lake temperatures and light levels increased with distance from the ice sheet, along with dissolved organic C (DOC) concentrations and total nitrogen:total phosphorus (TN:TP) ratios. Diatom production measured as biogenic silica accumulation rates (BSiARs) and diatom contribution to microbial communities declined across these gradients, while total production estimated using C accumulation rates and δ13C increased, indicating that autochthonous production and C burial are controlled by microbial competition and competitive displacement across physiochemical gradients in the region. Diatom production was generally low across lakes prior to the 1800’s AD but has risen 1.5–3× above background levels starting between 1750 and 1880 AD. These increases predate contemporary regional warming by 115–250 years, and temperature stimulation of primary production was inconsistent with paleorecords for ~ 90% of the last millennium. Instead, primary production appeared to be more strongly related to N and P availability, which differs considerably across the region due to lake landscape position, glacial activity and degree of atmospheric nutrient deposition. These results suggest that biological responses to enhanced nutrient supply could serve as important negative feedbacks to global change.

ABSTRACT Global change is reshaping the physical environment and altering nutrient dynamics across the Arctic. These changes can affect the structure and function of biological communities and influence important climate-related feedbacks (for example, carbon (C) sequestration) in biogeochemical processing hot spots such as lakes. To understand how these ecosystems will respond in the future, this study examined recent (< 10 y) and long-term (1000 y) shifts in autotrophic production across paraglacial environmental gradients in SW Greenland. Contemporary lake temperatures and light levels increased with distance from the ice sheet, along with dissolved organic C (DOC) concentrations and total nitrogen:total phosphorus (TN:TP) ratios. Diatom production measured as biogenic silica accumulation rates (BSiARs) and diatom contribution to microbial communities declined across these gradients, while total production estimated using C accumulation rates and d 13 C increased, indicating that autochthonous production and C burial are controlled by microbial competition and competitive displacement across physiochemical gradients in the region. Diatom production was generally low across lakes prior to the 1800's AD but has risen 1.5-39 above background levels starting between 1750 and 1880 AD. These increases predate contemporary regional warming by 115-250 years, and temperature stimulation of primary production was inconsistent with paleorecords for 90% of the last millennium. Instead, primary production appeared to be more strongly related to N and P availability, which differs considerably across the region due to lake landscape position, glacial activity and degree of atmospheric nutrient deposition. These results suggest that biological responses to enhanced nutrient supply could serve as important negative feedbacks to global change.

HIGHLIGHTS INTRODUCTION
Lakes are prominent features of Arctic landscapes (Paltan and others 2015). Here, they act as hot spots for nutrient cycling and carbon (C) burial that are highly sensitive to environmental forcing (Anderson and others 2009; Dean and others 2016). Lake ecosystem responses to global change drivers therefore represent key feedbacks to local-and global-scale biogeochemical dynamics; however, their relative importance is poorly understood due to an absence of long-term monitoring data. To fill this void, arctic lake sediment records have provided invaluable information on background environmental conditions and recent biological responses to anthropogenic change (Pienitz and others 2004). Much of this work has focused on shifts in algal (that is, typically diatom) community structure and production allegedly resulting from increased ice-free periods of lakes associated with high-latitude warming (Smol and others 2005;Burpee and Saros 2020). Although this mechanism undoubtedly plays a role in the rapid contemporary changes observed in many lakes, warming trends are highly asynchronous across the Arctic (Serreze and others 2000;Comiso and Hall 2014), and algal community responses can differ considerably at regional and local scales (Anderson 2000;Fortin and Gajewski 2009). Moreover, changes in energy inputs and mass loading into these ecosystems are occurring that could also drive shifts in algal communities and confound overly simplistic predictions of the temperature paradigm (Leavitt and others 2009;Saros and Anderson 2015).
In addition to temperature, primary production is fundamentally tied to light and nutrient supply in aquatic ecosystems. Ice and snow cover constrain algal growth in arctic lakes for most of the year by altering light quality/quantity and limiting external nutrient delivery (Roulet and Adams 1986;Weyhenmeyer and others 2013). Soon after ice-off however, intense solar radiation leads to rapid stratification and establishes lake mixing depths, which in turn control light regimes and mass/en-ergy flux, governing rates of algal photosynthesis (Diehl and others 2002;Saros and others 2016). Once harvested, light energy is used to fuel cellular metabolism and biomolecule synthesis (that is, nucleic acids and proteins) from essential elements such as C, nitrogen (N) and phosphorus (P) (Elser and others 1996). Narrow stoichiometric requirements for optimal growth exist for these elements (N:P 16 ± 3 mol), and because of this microbial N and P uptake/recycling is hypothesized to regulate their global-scale availability (Loladze and Elser 2011). These physiological constraints also mean that primary production is frequently limited by N and P (Elser and others 2007) and that changes in the supply of one or both elements can stimulate algal growth in oligotrophic arctic lakes (Levine and Whalen 2001;Hogan and others 2014).
Ice sheets and glaciers control physiochemical and ecological dynamics in paraglacial environments. Subglacial erosional and weathering processes of the Greenland ice sheet release enormous quantities of elements (for example, N at 30-40 Gg y -1 ; Wadham and others 2016 and P at 11 Gg y -1 ; Hawkings and others 2016) into surrounding ecosystems. Ice sheets also exert considerable influence on their surrounding environments by shaping local climate and topography. Legacies of Holocene deglaciation, beginning between 11 and 10 cal ky BP in areas of SW Greenland (Ten Brink and Weidick 1974), have established strong geochemical gradients from ice sheet margins to the coast that continue to shape biological structure and function in these areas (Anderson and others 2001;Deuerling and others 2019). More recently, regional glaciers generally retreated with higher temperatures prior to 1000 cal y BP followed by advances corresponding with cooling into the Little Ice Age (LIA) ending between the mid to late 1800's AD (Weidick 1984;Dahl-Jensen and others 1998;Carrivick and others 2018). Since then, they have undergone asynchronous minor advance/retreat cycles, the timings of which vary considerably at local and regional scales (Weidick 1984;Forman and others 2007). Beginning as early as the 1980's and as late as the 2000's AD, continuous retreats however have largely coincided with temperature increases and precipitation changes in the area (Forman and others 2007;Mouginot and others 2019). Integrated changes in glacial, terrestrial and aquatic ecosystems during this period have led to nonlinear shifts in both major environmental drivers and terrestrial and aquatic biological communities (Saros and others 2019).
The complexities of modern global change combined with legacy and contemporary effects of deglaciation impose substantial challenges for determining the relative influence of individual physiochemical drivers on lake ecosystem function. To address this problem, contemporary limnological and paleolimnological records were examined from eight lakes located along a 40 km transect originating near the present margin of the Greenland ice sheet. A replicated lake approach was used to encompass regional variability and determine (1) how physical and geochemical processes in paraglacial lakes differ across the landscape according to their proximity to the ice sheet. Spatiotemporal differences in diatom and microbial production were explored to (2) determine community and ecosystem-level responses to environmental change in the region over decadal and centennial time scales.

Study Area
Study lakes are small (< 25 ha) glacially scoured basins located in the interior of the Kangerlussuaq region between 2 and 41 km of the Greenland ice sheet (Figure 1). Freshwater lakes here are generally oligotrophic, chemically dilute and dimictic (Anderson and others 2001;Saros and others 2016), and under-ice light availability and primary production are influenced by snow depth (Whiteford and others 2016). Lakes stratify rapidly after ice-off, with mixing depths controlled predominately by light attenuation, although most lakes receive sufficient light to support phytoplankton production during ice-free periods (Saros and others 2016;Malik and others 2017). As there are no inflowing streams and groundwater inputs are negligible due to the presence of continuous permafrost, study lakes are largely hydrologically isolated from their surrounding catchments (Anderson and others 2001;Johansson and others 2015) and have been disconnected from the ice sheet for the last 7 ka (Carlson and others 2014). Annual precipitation values are low (< 250 mm y -1 ; Mernild and others 2015), and there is considerable sublimation of snowfall in this region resulting in minimal overland flows into the lakes (Anderson and others 2018). Due to hydrologic isolation, major biogeochemical inputs are thought to be predominantly derived from atmospheric fixation (C and N) and/or wet and dry deposition of N and P, respectively (Anderson and others 2017; Anderson and others 2018; Saros and others 2019).
Contemporary annual air temperatures are 2-3°C cooler, and monthly precipitation values are > 10 mm higher on average at higher elevation sites nearer the ice sheet (Table 1) compared to lowland sites at the head of the fjord (Whiteford and others 2016; Fowler and others 2020). Without detailed, site-specific historical temperature data, we assume that these general temperature gradients have held over the last millennium in our study area and that lakes have responded relatively synchronously to regional climate changes. Although these changes must be inferred from for most of the study period (1000-70 cal y BP; as described above), direct air temperature measurements began in Nuuk ( 330 km S of our study area) and Ilulissat ( 230 km N) during the 1880's AD at the end of the LIA. Regional temperatures did not begin rising however until between 1905 and 1915 AD (Box 2002), corresponding with local glacial advances up to 1920 AD. From this time, temperatures rose and local glaciers retreated into the 1940's (Knight and others 2000;Box 2002), which is when long-term temperature records become available for Kangerlussuaq. Local glaciers continued retreating until the late 1960's and early 1970's despite cooling trends in the area beginning in the early 1940's (Van Tatenhove 1996; Knight and others 2000). Following this delay, local glaciers advanced with cooling temperatures into the mid to late 1990's (Knight and others 2000; Saros and others 2019) with both trends reversing rapidly with regional warming into the present.

Sample Collection, Processing and Analyses
Three pairs of lakes were sampled for this study to provide regional replication across the landscape. Sites were chosen based upon their landscape position ( Figure 1) with SS906 and SS901 referred to as ''Ice Sheet'' lakes, and lakes 17b and 17c found in the middle of the study region denoted as ''Central.'' Lakes SS1590 and SS85 are located furthest away from the ice sheet near the head of the fjord. We refer to these sites as ''Kellyville'' lakes, using a common local name, while noting that these lakes have previously been referred to as inland or central sites in studies extending out to the coast (Whiteford and others 2016; Osburn and others 2019). Data from two previously studied lakes (SS1341 and SS904; Anderson and others 2018) were also included to extend the spatial extent of contemporary analyses.
Sampling and processing of water physiochemical parameters followed protocols used in previous limnological surveys of the area (Anderson and others 2001; Whiteford and others 2016; Curtis and others 2018; see Supp. Methods S1 for exhaustive details). Water nutrient fractions including: total P (TP), dissolved inorganic N (DIN) and silicate (SiO 2 ) were analyzed in a temporary laboratory established at the Kangerlussuaq International Sciences Support (KISS) station using common standardized methods (APHA 2017). Water total N samples were frozen at -20°C and shipped to the CRASR laboratory at Baylor University (USA) for analysis. Dissolved organic C samples were acidified and brought back to Loughborough University (UK) for analysis. Water physiochemistry and landscape data from the present study were combined with those from the previous work to synthesize a comprehensive suite of variables that could influence autochthonous production in the region (see Table 1). Sediment cores were taken from the deepest part of each lake using a HON-Kajak gravity corer (Renberg 1991). Cores were extruded on site at a resolution of 0.25 to 0.5 cm intervals. Sediments were frozen at -20°C and brought back from Greenland for subsequent analyses. Bulk density was determined for each sediment interval by drying a known volume of sediment overnight at 100°C; afterward, organic matter content was determined according to Dean (1974). Remaining lake sediments were then lyophilized and sent for 210 Pb-dating at the St. Croix Watershed Research Station using an EG&G nuclear ultra-low background gamma spectrometer (Ortec; TN, USA). Sediment dating methods of the other two previously collected cores are comparable and discussed elsewhere (Anderson and others 2018).
Chronologies for 210 Pb-dated parts of all cores were derived using the constant rate of supply model (Appleby 2001). For the six newly collected paired lake cores, 14 C dating was conducted using accelerator mass spectrometry at the Natural Environment Research Council (NERC) Radiocarbon facility (East Kilbride, UK). Two bulk sediment samples were submitted for 14 C dating from four out of six lakes (SS1590, 17b, 17c and SS906), with one sample being submitted for previously dated lakes (SS85; Olsen and others 2013 and SS901; Whiteford and others unpublished). One additional date for each lake was determined through crosscorrelation of LOI profiles from previously dated cores (Supp. Figure S2). Radiocarbon dates were calibrated and combined with 210 Pb-dates (see Supp. Tables S3-8 for dating information) to produce a single age-depth model for each lake using the R ''clam'' package v.2.3.2 (Blaauw 2010).
Sediment elemental and isotope analyses were performed using isotope ratio mass spectrometry and inductively coupled plasma mass spectrometry according to Osburn and others (2019) and Watts and others (2019), respectively. Sediment BSi concentrations were measured on alternating sediment layers using standard 1% sodium bicarbonate extractions and silicate described in Conley and Schelske (2001). Diatom production, reported as BSi accumulation rates (BSiAR), was calculated by multiplying sediment %BSi by dry mass accumulation rates for discrete sediment intervals. To standardize BSiARs for cross-lake comparisons, dry mass ARs were focus-corrected by calculating ratios of regional atmospheric 210 Pb flux and mean 210 Pb flux in individual sediment cores using previously reported mean atmospheric flux rates for the region (Anderson and others 2019). It is acknowledged that BSi is a measurement of total siliceous algal biomass including both diatoms and chrysophytes. However, since it is impossible to distinguish between the two taxa using BSi alone, we refer to BSiAR as a measurement of diatom production for simplicity and consider the distinct contemporary ecological differences between the two groups in our landscape BSiAR interpretations.

Statistical Analyses
Visual inspection of temporal changes in BSiARs and BSi:C ratios in paired lake cores revealed nonlinear dynamics, and as a result, these patterns were assessed using general additive models (GAMs) fit using restricted maximum likelihood methods using the ''mgcv'' v.1.8-31 package (Wood 2011). Thin plate regression splines and ktests reported in the gam.check function were used to optimize knot selection for each lake (n = 9), and temporal autocorrelation was accounted for using a continuous-time first-order autoregressive process, according to Simpson (2018). To detect periods of significant changes in response variables, simultaneous confidence intervals were calculated using the ''gratia'' package (Simpson 2018). Significant changes in diatom production were estimated from the first derivative of fitted trends using the method of finite differences.
For spatial trend analyses, periods before and after BSiAR increases in each lake were identified from GAMs, and BSiAR enrichment ratios were calculated as the ratio of pre:post-BSiAR estimates. One-way ANOVAs on log-transformed values were used to determine significant differences between lake groups identified from enrichment ratios. Significant differences among sites were then compared using Bonferroni corrected Tukey's post hoc tests. To examine correlations between landscape parameters, contemporary water column physiochemical variables and surface sediment BSiAR, Pearson correlations were used for all monotonic relationships. A unimodal relationship between mean depth and BSiAR was modeled using a 3-parameter Gaussian function that was selected from a suite of potential nonlinear functions using Akaike information criterion.
In addition to contemporary analyses, centennial-scale relationships between sediment molar N:P ratios, autotrophic and microbial production (BSiAR and d 13 C), carbon burial (CAR) and diatom relative biomass contribution to microbial assemblages (%BSi:C) across paired lakes were examined using Pearson correlations after linear-detrending of each parameter and log or pseudo-log transformation (d 13 C only). Inspection of residual plots revealed 1-3 outlier points for each plot, defined as greater than 2 SD from mean lake values, so these points were removed, and correlations reran. Similar methods were used to examine correlations between sediment stoichiometry and BSiAR within each lake. However, BSiAR outliers, which were found in 17c and SS901, were retained for these analyses out of necessity due to the smaller sample sizes of individual lakes. Correlations in these two lakes were determined using nonparametric Kendall tau tests.
As global and regional algal production is often limited or co-limited by N and P (Burpee and others 2016; Elser and others 2007; Whiteford and others 2016), changes in algal biomass associated with lake N and P availability have traditionally been used as limnological indicators of short-term nutrient limitation. Since BSiAR serves as a proxy for diatom, and more generally autotrophic production in lakes (Conley and Schelske 2001), mass balance principles can be used to extend these inferences to longer timescales by examining correlations between BSiAR and sediment nutrients, provided that certain assumptions regarding deposition and preservation of these signals within lake sediments are met (see Interlandi and others 2003; Supp. Text S9).

RESULTS
Spatial patterns in contemporary (< 10 years) BSiARs were related to lake landscape position and associated regional physiochemical gradients (Table 1). Sediment BSiARs decreased with distance from the ice sheet ( Figure 2; r = -0.76, p = 0.030) and were positively correlated with altitude (r = 0.93, p < 0.001) with the highest accumulation rates ‡ 9.5 g m -2 y -1 recorded at the two highest elevation lakes located nearest to the ice sheet. Diatom production was nonlinearly related to mean depth (r = 0.83, p = 0.056), although this trend could be considered marginally significant. Lake water conductivity, DOC and TN increased with distance from the ice sheet, and these patterns were associated with decreased BSiAR (p < 0.05). There was a positive trend between water TP and BSiAR (r = 0.62, p = 0.100), and an associated negative relationship between TN:TP ratios and BSiAR (r = -0.79, p = 0.021) identified strong linkages between proportional water N and P availability and diatom production across the region. No other significant relationships with physiochemical variables were detected (water temperature, lake area, max depth, pH, alkalinity, DIN and SiO 2 ; p > 0.14).
Landscape Controls on Arctic Lake Production Looking at historical patterns in paired lakes, focus-corrected sediment BSiAR values ranged from 0.1 to 6.5 g m -2 y -1 from 1000 to -67 cal y BP (Figure 3). Background BSi accumulation rates were relatively low (0.6 ± 0.3 g m -2 y -1 ) and similar across lakes throughout most of the sediment record until conspicuous increases (1.5-fold to threefold) began 200 cal y BP (Figures 3 and  4). After this time, lakes could generally be separated by landscape characteristics with lower BSiARs (0.8 ± 0.3 g m -2 y -1 ) in lower elevation lakes located closer to the fjord (hereafter referred to as ''low production'' lakes; SS85, SS1590 and 17b) compared to ''higher production'' lakes found at higher elevations closer to the ice sheet (17c, SS901 and SS906; 2.1 ± 1.4 g m -2 y -1 ; Figure 4). These groupings reinforce contemporary produc-tion patterns and were also apparent at the centennial scale as BSiAR in low production lakes generally showed unimodal changes from 1000 to 200 cal y BP followed by exponential BSiAR increases, whereas more productive lake values remained stable during this period before increasing exponentially after 200 cal y BP ( Figure 3).
As with contemporary water column measurements, proxies of diatom and microbial production were related to N:P ratios, which increased with distance from the ice sheet. Sediment BSiAR and d 13 C decreased while CAR increased significantly across sediment N:P gradients (p £ 0.001; Figure 5), indicating reduced diatom production and greater C burial in lakes located further from the ice sheet. Declines in sediment %BSi:C ratios across N:P gradients (r = -0.72, p < 0.001) along with higher sediment CAR (r = -0.29, p < 0.001) suggest that these changes are due to reduced contribution of diatom biomass to the sediment C pool. Centennial-scale relationships between sediment nutrients and BSiAR provide evidence of differential nutrient limitation among study lakes. Kellyville lakes located furthest from the ice sheet, SS85 and SS1590, appeared to be predominantly P limited (Table 2; Supp. Figure S10). Lake 17b in the center of the study region showed signs of N and P co-limitation, but no evidence of elemental limitation in 17c was found. The lakes closest to the ice sheet, SS901 and SS906, both appeared to be colimited with no changes in sediment %N, P or BSiAR observed in their sediments before 200 cal y BP (Supp. Figure S11).

DISCUSSION
Strong climate and geochemical gradients exist in paraglacial environments that influence ecosystem responses to global change (Matthews 1992;Anderson 2007;Engstrom and others 2000). To predict future changes in these sensitive areas, a better understanding of past and present biological responses to environmental variability across these landscapes is needed. This study addresses this need by examining contemporary patterns and centennial-scale changes in autotrophic production and interpreting these results through insights from extensive limnological research conducted in the Kangerlussuaq region of SW Greenland.

Physiochemical Drivers of Contemporary Spatial Variation in Diatom Production
Diatom production was correlated with many physical and chemical variables related to lake landscape position (Figure 2). Contemporary BSiARs decreased with distance from the ice sheet and were positively correlated with altitude. Therefore, unlike most widely held assumptions of strong relationships between diatom production and temperature (Douglas and Smol 1999), BSiARs were reduced at warmer (2-3°C on average) lowaltitude sites and showed no correlation with epilimnion water temperatures. Although seasonal air temperatures are correlated with ice-off dates and phytoplankton biomass across the region (Whiteford and others 2016), counter gradient responses in diatom production suggest that temperature is most likely an indirect factor colinearly associated with other key limnological controls on BSiARs. In support of this hypothesis, diatom growth rates are unrelated to independent effects of water temperature in the area, and only with adequate light and nutrient supply does temperature appear to play a role (Saros and others 2014;Malik and Saros 2016). Nevertheless, given the rapid nonlinear changes in temperature and a suite of physiochemical changes in the region over the last 30 years (Saros and others 2019), temperature is likely to become a more important direct or indirect factor related to regional primary production in the future as has been reported in other high-latitude and alpine lakes (Moser and others 2019; Burpee and Saros 2020).
Diatom production was also negatively correlated with DOC concentration. As autotrophs, DOC availability does not directly affect diatom production but instead could influence BSiAR by controlling the light climate. Waters with elevated DOC and colored dissolved organic matter (DOM) concentrations absorb photosynthetically active radiation (PAR) in arctic lakes resulting in shallower mixed layers and higher average epilimnion light levels (Saros and others 2016). Given adequate nutrient supplies, some pelagic diatom species (for example, D. stelligera) and chrysophytes can thrive under low-light conditions, and experimental evidence along with regional phytoplankton surveys suggests that diatoms could be competitively displaced by other phytoplankton taxa such as dinoflagellates and cyanobacteria under high light/DOC conditions found in many Kellyville lakes (Malik and Saros 2016; Northington and others 2019). Lake DOC and light availability are positively correlated with water temperature across the study area, and these variables interact to control lake mixing depths in Kangerlussuaq lakes (Northington and others 2019; Saros and others 2016). However, because the effects of light on mixing depths are considerably greater than those of temperature (Saros and others 2016), light seems to represent the predominant physical control on diatom production in the region.
Light availability can also interact with lake morphology to control diatom production rates. Mean depth and depth:PAR ratios influence diatom production by affecting the light climate of pelagic areas and determining the amount of benthic area available for diatom colonization (Vadeboncoeur and others 2008;Whalen and others 2008;Malik and others 2017). In lakes with shallower mean depths located furthest from the ice sheet, low BSiAR values can be attributed to competitive displacement of pelagic diatoms in high light environments as outlined above (Northington and others 2019). Further, the presence of extensive cyanobacterial mats in the shallowest study lake SS85 (Anderson personal observation) suggests that similar community-level processes might also be operating in benthic habitats. The two deepest lakes on average, 17b and 17c, had similarly low contemporary BSiARs, which could be tied to lightlimitation of benthic diatoms across large areas of these basins. The highest BSiARs were found in ice sheet lakes where relatively deeper mixed depths and lower average light availability may favor diatom dominance in pelagic algal communities (Malik and others 2017; Northington and others 2019; Saros and others 2016), along with the potential for moderate benthic production at intermediate mean depths. These patterns reinforce the importance of benthic environments for autotrophic production in Arctic lakes (Whalen and others 2008); however, while light structures phytoplankton assemblages across the region, most Kangerlussuaq lakes receive sufficient PAR at mean depth to support phytoplankton growth (Malik and others 2017). This suggests that light availability primarily controls regional diatom production by influencing microbial competition and community composition rather than through physiological light-limitation per se.

Nutrient Controls on Contemporary Regional Diatom Production
Lake nutrient availability also regulates algal community composition patterns across the Kangerlussuaq region (Hogan and others 2014;Malik and others 2017;Northington and others 2019). Externally supplied DIN is primarily delivered through precipitation (Anderson and others 2018), which is slightly higher at the ice sheet compared to Kellyville (Fowler and others 2020). However, DIN is only a minor component of TN pools and does not differ significantly across the study gradient (Table 1; Anderson and others 2018; Curtis and others 2018). Lake TN is related to DOC (and conductivity) due to strong effective precipitationevaporation balance gradients across the region (Anderson and others 2001), suggesting that N-rich DOM is the predominant source of N in Kangerlussuaq lakes (Anderson and Stedmon 2007;Burpee and others 2016). As much of the regional DOM pool is composed of refractory material (Anderson and Stedmon 2007), a portion of this organic N is unavailable or metabolically costly to acquire for diatoms, which could explain their competitive displacement by N-fixing cyanobacteria and mixotrophic algae (that is, dinoflagellates) in DOM-rich Kellyville lakes (Northington and others 2019). Conductivity and DOM gradients could also explain changes in silicious algal taxonomic composition across the region, as mixotrophic chrysophytes can leverage organic C sources and fair better than diatoms in more highly conductive waters (Northington and others 2019). In contrast to chrysophytes, diatom proportional biomass is much greater in lakes located closer to the ice sheet and sandurs (Figures 1 and 2; Supp. Figure 12), which supply considerable amounts of P-rich sediments to the surrounding area (Hawkings and others 2016; Storms and others 2012). Katabatic winds blowing away from the ice sheet produce relatively frequent and substantial dust storms in Kangerlussuaq (Bullard and Mockford 2018), but these winds tend to weaken 10 km away from the ice sheet (Heinemann 1999). This suggests that decreasing dust transport to areas further away from the ice sheet could partially explain trends of lower lake TP and P-limited microbial growth in Kellyville lakes (Brutemark and others 2006;Burpee and others 2016). Shifts from pelagic to benthic production, or ''benthification,'' that have occurred in Kangerlussuaq since the mid-1990's AD are likely to exacerbate P deficiencies in these lakes through biological sequestration in littoral areas (Saros and others 2019).
Although N and P availability regulates algal community structure in the region, resource stoichiometry determines the nature of diatom nutrient limitation (Tilman 1977). Lake DIN:TP (and light:nutrient) ratios shape relative diatom species abundances across Kangerlussuaq lakes (Malik and Saros 2016;Malik and others 2017), but these ratios are unsuitable for predicting microbial nutrient limitation and community composition in the area (Burpee and others 2016; Northington and others 2019). This is likely due to the greater proportion of DON to DIN in these lakes (Anderson and others 2018), meaning that TN:TP ratios could be a more accurate predictor of both diatom relative biomass and potentially their nutritional status. Contemporary BSiARs were the highest in ice sheet lakes with TN:TP ratios below 50 with a distinct drop off in diatom biomass observed in the next two most productive lakes with TN:TP of 75 and 100 (Figure 2). Year-round microbial P-limitation has been reported in pelagic zones of many of the higher TN:TP Kellyville lakes (Brutemark and others 2006;Burpee and others 2016), but more data are needed to precisely determine N:P threshold elemental ra-tios of regional diatom nutrient limitation. Regardless, this evidence indicates that contrary to the widely held paradigm of the predominance of N-limitation at higher latitudes (Vitousek and Howarth 1991), P-limitation and N and P co-limitation are important constraints on autotrophic production in the Kangerlussuaq region (Burpee and others 2016; Hogan and others 2014) and potentially across the Arctic (Levine and Whalen 2001;Symons and others 2012). Extending these contemporary insights back over greater time periods can yield a more detailed understanding of how landscape factors and climate control arctic diatom production by mediating stoichiometric nutrient availability.

Landscape, Climate and Anthropogenic Controls on Centennial-Scale Lake Production
Longer-term spatial patterns in diatom and microbial community production across lake light/nutrient gradients highlight the role of autochthonous production in shaping regional biogeochemical cycling. Diatom production was uniformly low across the study area until exponential increases began in each lake between 200 and 70 cal y BP (Figure 3). From the start of this period, lakes can be separated into lower production lakes located closer to the head of the fjord and higher BSiAR lakes located nearer to the ice sheet (Figure 4b). Differences in sediment BSi:C ratios across these lakes confirm the importance of competitive displacement in structuring regional algal assemblages, as diatom relative biomass decreases proportionally with increasing lake N:P ratios ( Figure 5). Lower sediment d 13 C and BSiAR's contrasted with higher CAR values across these gradients would seem to indicate that diatom and microbial community production are negatively related across the study area. However, this interpretation must be tempered by observations of strong collinearity in physiochemical gradients across the region that can influence lake C dynamics, greater C inputs from vascular plants in Kellyville lakes and high minerogenic inputs of glacial dust closer to the ice sheet (Anderson and others 2012, 2019; Osburn and others 2019). As BSiAR and CAR are commonly used as independent proxies of whole algal community production, understanding discrepancies in these rates and quantifying contributions of both diatom and other sensitive algal taxa (for example, cyanobacteria and chlorophytes) to C burial across landscape gradients is essential and requires detailed studies linking Landscape Controls on Arctic Lake Production Figure 3. Temporal changes in focus-corrected biogenic silica accumulation rates (BSiAR). Discrete measurements from lake pair sediment cores (n = 6) are depicted as black circles. Sediment BSiARs (solid lines) were modeled using general additive models with dashed lines indicating 95% confidence intervals. All model fits were statistically significant (p < 0.001) with adjusted-R 2 values ranging between 0.56 and 0.76. Changes in BSi accumulation rates were determined using confidence interval derivatives and are indicated with vertical dashed lines (and shown in red in the online version). Temporal BSiAR thresholds (cal y BP) for each lake are: A SS85, 78; B SS1590, 128; C 17b, 92; D 17c, 135; E SS901, 144 to -37; and F SS906, 198. *Note y-axes are shown on different scales to highlight temporal variation within lake pairs for Kellyville (A, B), Central (C, D) and ice sheet (E, F) lakes where distance from ice sheet decreases from (A--F).
water column metabolism measurements to community composition data.
Nevertheless, BSiAR and CAR are strongly correlated within individual lakes in this region (r = 0.78-0.97), meaning that a replicated regional approach can provide valuable information on environmentally controls on algal production and C burial. Prior to 1000 cal y BP, local glaciers receded with higher temperatures most likely during the Medieval Climate Anomaly (Forman and others 2007;Carrivick and others 2018). Kangerlussuaq glacial dynamics are poorly constrained from 1000 to 500 cal y BP, but mass balance modeling in SW Greenland and ice-core temperature reconstructions suggest glacial advances with neoglacial cooling during this period ( Figure 6A; Dahl-Jensen and others 1998;Weidick 1984). Precipitation decreased considerably from 700 cal y BP (Aebly and Fritz 2009), which coincided with increased winds and/or glacial sediment availability, resulting in enhanced regional dust deposition peaking between 550 and 450 cal y BP ( Figure 6B; Anderson and others 2012;Willemse and others 2003). Increased BSiAR in low production lakes and stable rates in most likely N-limited high production lakes from 1000 to 400 cal y BP suggest that P loading from glacial dust was the primary factor controlling diatom production in lakes closer to the fjord during this time. As temperature was either negatively or neutrally related to diatom production across all lakes, thermodynamics dictate that any temperature effects on BSiAR were likely indirect and related to wind/glacial activity since increased ice cover with lower temperatures would have reduced light inputs and slowed diatom metabolism during this period.
Differential responses among lower production lakes underscore the importance of accounting for landscape heterogeneity in paleolimnological reconstructions (Anderson 2014). After initial cooling, climate began to warm following an initial LIA minimum 400 cal y BP (Dahl-Jensen and others 1998), and dust deposition rates decreased until after 200 cal y BP resulting in unimodal patterns in BSiAR from 1000 to 200 cal y BP in SS1590 and 17b (Figures 3 and 6). In contrast, BSiAR in SS85, which exhibited a more muted peak in dust deposition during this time (Anderson unpublished data), declined slightly. The distance of SS85 from the ice sheet likely resulted in reduced dust inputs and sustained P-limitation, partially explaining lower diatom production rates in SS85 compared to other lakes. Differences in diatom contributions to microbial production could also have played a role as BSiARs and BSi:C ratios declined synchronously Figure 4. Spatial differences in biogenic silica accumulation rates (BSiARs) over the last millennium. Median BSiARs are shown for each lake (n = 6), along with 25th and 75th percentiles (boxes), 10th and 90th percentiles (whiskers) and outlier data points (black circles). Baseline lake BSiARs are shown between A 1000 BP and the onset of BSi increase within each lake, along with B Post-BSi increase values for the remainder of the time series. Enrichment ratios C highlight increases in BSiAR over baseline values. Differences in BSiAR among individual lakes were determined through one-way ANOVA's followed by Bonferroni corrected Tukey's post hoc comparisons and are denoted by letters. Significant BSiAR differences were found between low and high production lakes both before and after BSiAR increases (p < 0.001). Significant differences between pre-BSi increases were also found between Central lakes (17b and 17c), whereas Kellyville (SS85 & SS1590) and ice sheet (SS901 & SS906) lakes different post-enrichment. in SS85 while showing unimodal changes in SS1590 between 1000 and 200 cal y BP (Figure 3 and Supp. Figure S12). Because opposite patterns in diatom relative biomass would be expected with changes in temperature, nutrient supply rather than physical constraints seems to have had a greater influence on diatom production in lower production lakes over most ( 90%) of the last millennium. Further, BSiARs across all lakes were much lower from 1000 to 200 cal y BP in comparison with more recent exponential increases across the region, suggesting an abrupt change in nutrient supply and ecosystem functioning after this time.
Climate effects on glacier dynamics have noticeably altered the generation/transport of P-rich sediments in the region over the last 250 years. Local glaciers reached their maximum late Holocene extents of up to 1-2 km for the Isunguata Sermia glacier and 100-200 m at Russell Glacier ahead of present-day ice margins toward the end of the LIA (245-70 cal y BP; Carrivick and others 2018;Forman and others 2007). From this point onward, multiple minor glacial advances and retreats have occurred, which have increased the quantity of sediments released to downstream sandurs (Weidick 1984;Knight and others 2000;Forman and others 2007) and potentially changed sediment quality (that is, nutrient content/lability) by producing newer and less weathered materials (Anderson 2007;Taylor and Blum 1995). As dust transport in Kangerlussuaq is primarily limited by Figure 5. Regional relationships between sediment nutrients and autotrophic production over the last millennium in paired study lakes (n = 6). Significant correlations (p £ 0.01) are shown for: A sediment nitrogen:phosphorus (N:P) ratios and focus-corrected biogenic silica accumulation rates (BSiARs), r = -0.44; B sediment N:P ratios and d13C, r = -0.27; C sediment N:P and BSi:carbon (C) ratios, r = -0.72 and D sediment carbon accumulation rates (CAR) and BSi:carbon (C) ratios, r = -0.29. Outlier data points (> 2 SD of lake means) are shown in white but were omitted from correlations. Correlations between sediment N:P and CAR, r = 0.38 and CAR and d13C, r = -0.56, were also significant (p < 0.001; data not shown). sediment availability (Bullard and Austin 2011), synchronous increases in P-rich dust inputs and lake BSiAR following the LIA highlight tight linkages between lake primary producers and the surrounding paraglacial environment.
In addition to P, N appears to become an important element for diatom production closer to the ice sheet. Following an initial production peak 600 cal y BP, BSiAR in 17b declined and did not begin increasing again until the second LIA minimum ( 1850 AD; Dahl-Jensen and others 1998), only surpassing its previous millennial maximum 1950 AD during a period of regional cooling (Anderson and others 2019). This timing corresponds with a noticeable increase in sediment %N (Supp. Figure S10), coinciding with the onset (inferred from D 15 N values; Hastings and others 2009) and peak increase in anthropogenic N deposition in the region according to ice-core NO 3 records ( Figure 6; Hastings and others 2009;Wolff 2013). Therefore, long-range N transport combined with high dust P inputs could explain more recent increases in diatom production in this co-limited lake. In contrast, there was no conclusive evidence of nutrient limitation in 17c where sediment C and N were unrelated to P, indicating that dilution effects of minerogenic dust could have obscured signs of nutrient limitation in this lake. It is also possible that other factors such as light exert greater control over diatom production in 17c due to its high mean depth and colored DOM concentrations (Anderson unpublished data), underscoring the importance of using coupled energy-mass balance perspectives when inferring drivers of environmental change from paleorecords (Leavitt and others 2009).
Production increases within ice sheet lakes began 50-100 years prior to the onset of anthropogenic N deposition (Figure 6), suggesting that locally derived N must also be important for autotrophic production in these lakes. Lake catchment inputs are unlikely to have been an important factor as low regional temperatures and precipitation rates would have limited groundwater inputs and restricted overland flow and nutrient delivery at this time. However, a potentially overlooked regional N source-the surrounding ice sheet-contains enormous amounts of DIN (40 tg; Wolff 2013). Local glaciers would have been close to their LIA maximum when %N and BSiAR began to increase in ice sheet lake sediments ( 1800 AD). Volatilization or photolysis of ice sheet DIN (Zatko and others 2016) combined with stronger katabatic winds and a lower boundary layer under colder temperatures (Beare and others 2006) could have provided enough N to partially alleviate N-limitation in these lakes. Similarly, photolysis of greater amounts of snow NO 3 during the LIA could have been another important local supply of N as has been observed in both SW Greenland and Antarctic regions (Curtis and others 2018;Wolff 2013).
Although the relevance of external N-sources remains to be determined, there is direct observational evidence that N fixation could explain diatom production increases in ice sheet lakes prior to the onset of long-range N deposition. Increased cyanobacterial biomass has been documented from 200 cal y BP in SS901 (Whiteford and others unpublished). Here, high N fixation/mineralization rates could stimulate diatom production, especially within benthic areas, which can account for around half of autotrophic production in lakes (Vadeboncoeur and others 2002;Whalen and others 2008). An even more direct connection is found in the diatom Epithemia sp., which hosts an N-fixing symbiont and is found in high contemporary abundances in regional ice sheet lakes (McGowan and others 2018). Despite the potential role of localized N-inputs however, dramatic increases ice sheet lake BSiAR over the last two centuries are unlikely to have been supported by local N alone. Overall, these findings along with extensive aeolian activity documented in the re- gion demonstrate the importance of both longrange and localized N and P subsidies for supporting autotrophic production in this region and likely for many other parts of the Arctic (Anderson and others 2017). By combining insights from paleo-and contemporary limnology, this study outlines two common mechanisms-competition and nutrient limitation-that structure primary producer assemblages in SW Greenland. These factors are likely to influence algal ecology across high-latitude lakes, but results from this study and others show that, even within small geographic regions, lake ecosystems can respond quite differently to changes in climate and nutrient loading (Anderson 2000;Fortin and Gajewski 2009). Therefore, wider extrapolation requires detailed consideration of hierarchical global, regional and local-scale processes. If autotrophic communities in Kangerlussuaq lakes continue operating as they have over the last 250 years, their biogeochemical role in offsetting rises of global greenhouse gas concentrations through increased production and C burial (Anderson and others 2019; Anderson and others 2009), will only continue to become more important. However, rapid warming since the turn of the century has significantly altered local glacier mass balance, riverine discharge and dust events in Kangerlussuaq (Bullard and Mockford, 2018; Mouginot and others 2019; Saros and others 2019) coinciding with rapid and unpredictable changes to terrestrial and aquatic ecosystems (Saros and others 2019). This increased variability in energy and mass flux will undoubtedly have unforeseen effects on lake structure and function over the coming years, but detailed studies incorporating stoichiometric concepts could help parse through this complexity and better resolve the biogeochemical roles of lakes in a warmer Arctic. Figure 6. Changes in temperature, nutrient deposition and diatom production over the last millennium. A Harmonized ice sheet temperature trends are modeled with z-scores from GRIP and DYE 3 ice cores (Dahl-Jensen 1998) using a general additive model (GAM). B Focused-corrected minerogenic accumulation rates (FCMinARs), which serve as a proxy for regional dust and phosphorus deposition, are modeled separately in low production and high production lakes using GAMs. Ice sheet nitrate (NO 3 ) values are replotted from an ice core collected near the Greenland Summit (B16; Wolf 2013). Focus-corrected biogenic silica accumulation rates (FCBSiARs) are modeled for 3 C low and D high production lakes using GAMs where solid lines show production trends, dashed lines denote 95% confidence intervals and vertical dashed lines, and solid red lines indicate periods of significant change C 100 cal y BP D 164 cal y BP. Polar Field provided invaluable logistical support. Two anonymous reviewers supplied helpful feedback that improved the manuscript from its original version.

OP E N A C C ES S
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 h ttp://creativecommons.org/licenses/by/4.0/.