Changes in coupled carbon‒nitrogen dynamics in a tundra ecosystem predate post-1950 regional warming

Arctic ecosystems are changing in response to recent rapid warming, but the synergistic effects of other environmental drivers, such as moisture and atmospheric nitrogen (N) deposition, are difficult to discern due to limited monitoring records. Here we use geochemical analyses of 210Pb-dated lake-sediment cores from the North Slope of Alaska to show that changes in landscape nutrient dynamics started over 130 years ago. Lake carbon burial doubled between 1880 and the late-1990s, while current rates (~10 g C m−2 yr−1) represent about half the CO2 emission rate for tundra lakes. Lake C burial reflects increased aquatic production, stimulated initially by nutrients from terrestrial ecosystems due to late-19th century moisture-driven changes in soil microbial processes and, more recently, by atmospheric reactive N deposition. These results highlight the integrated response of Arctic carbon cycling to global environmental stressors and the degree to which C–N linkages were altered prior to post-1950 regional warming. Biological productivity and carbon dynamics in Arctic ecosystems started to change prior to human-induced warming of the region, according to an investigation of coupled carbon–nitrogen cycle dynamics using stable isotope analyses of lake sediments.

T errestrial carbon stores in the Arctic are substantial 1 and are mainly held in permafrost 2 ; as such they represent a potential significant feedback into the climate system as soils warm 3,4 . However, as many key processes (e.g., photosynthesis, soil microbial respiration) in arctic ecosystems are often nitrogen (N) limited 5 , future alterations to organic carbon (C) pools may be linked to changes in N fluxes 6,7 , including fixation, cycling, and transfer among adjacent ecosystems 8 . Recent warming at high latitudes has resulted in deepening of the seasonal active layer and permafrost thawing with an associated increase in mineralisation of organic matter and release of CO 2 and CH 4 to the atmosphere 4 . Such warming also alters soil microbial processes, which releases nutrients (especially N 9 ) that stimulate both terrestrial 10 and aquatic productivity 11 , with potentially important consequences for regional rates of C storage.
Future impacts of changing temperature and precipitation on nutrient cycling and ecosystem function are difficult to assess because of the inherent limitations of short-term experimental approaches and the paucity of decadal-scale monitoring records, particularly at high latitudes. Similarly, few studies have evaluated how the coupling between C and N cycles may have shifted in response to past climatic variability 12 . As lakes are a key component of arctic landscapes and are often net heterotrophic (net ecosystem respiration > primary production) 13,14 , they may be critical sites to quantify centennial-scale responses of C and N biogeochemistry to climate change. For example, atmospheric warming is expected to increase CO 2 efflux from lakes 15 ; however, associated changes in terrestrial C storage and N cycling as well as altered hydrological pathways 16 may also alter lateral transfer of terrestrial organic carbon as particulate (POC) and dissolved (DOC) to streams and lakes 17 , where it may be either stored in sediments 14 or released to the atmosphere following bacterial respiration or photochemical mineralisation 18,19 . The degree to which these processes offset each other are poorly constrained 17,20 , in part because of complex changes in mobilisation and interaction of C and N as permafrost thaws 21 .
The landscapes of the North Slope of Alaska are now undergoing pronounced permafrost thaw (deepening of the active layer), altered seasonality, and declining snow cover 22,23 . While atmospheric temperatures have increased terrestrial productivity, warming soils have enhanced soil microbial activity and respiration thereby elevating CO 2 and CH 4 release 4 . At the same time soil N dynamics 24 , shrub expansion (greening) 25 , thermokarst, and soil erosion are all changing in complex manners 22 . These terrestrial processes have cascading effects on lakes and wetlands 14 , which are important landscape components for both C processing and storage 18,26,27 .
Here, we use whole-basin organic C burial rates and stable isotope analyses of C and N (δ 13 C; δ 15 N) in 210 Pb-dated sediment cores (see Supplementary Fig. 1) to quantify how historic climate variation has altered coupled C and N cycles in five lakes located near the Toolik Lake Long-Term Ecological Research (LTER) Field Station in the Kuparuk River Basin, Alaska, USA 28 (Table 1). Data from these sediment records (see Supplementary Data #1) provide an integrated record of lake and catchment nutrient dynamics over the last~200 years, and show that a combination of changing moisture availability, climate warming, and reactive N deposition over the past 130 years has increased aquatic production and aquatic C burial regionally. Importantly, changes to lake and catchment nutrient cycling were initiated in the late-nineteenth century by a post-Little Ice Age climatic rebound 29 well before the rapid and more widely recognised post-1950s warming.

Results and discussion
Aquatic and sediment carbon dynamics. Sediment records from the five study lakes show that mean regional C burial rates slowly increased from~4 g C m −2 yr −1 in the nineteenth century, more than doubling by the 1990s (Fig. 1c). Both rates and trends over time are markedly similar across all sites ( Supplementary Fig. 2). The background burial rates (~4 g C m −2 yr −1 ) characterise the natural C sequestration rate for lakes situated in moist, acidic tussock tundra, while present-day rates (~10-12 g C m −2 yr −1 ) represent nearly 50% of total C emissions from freshwater lakes in this region (~24 g C m −2 yr −1 ) 18 . Most importantly, burial rates increased substantially during a period of 100 years.
There is also clear sedimentary evidence that the increased C burial rates arose mainly from enhanced in-lake production, driven by changes in lake catchment nutrient cycling and increased N influx. Present-day sediment C/N ratios across all sites are <15 and δ 13 C values are around -30‰, both of which are essentially algal signatures 30 and substantially lower than corresponding values for tussock peat or catchment runoff (Fig. 2a). In addition, the temporal trends in the sediment δ 13 C data are broadly similar ( Fig. 1c and Supplementary Fig. 2b) and suggest that aquatic primary producers used CO 2 respired from terrestrially fixed C rather than atmospheric CO 2 13 . This process is recognised as being important today 19 , but the δ 13 C profiles suggest that C dynamics started to change from~1880 (Breakpoint = 1877; see "Methods" and Supplementary Fig. 3 for summary of breakpoint analyses). Algal communities-as summarised by PCA axis-1 scores of fossil diatom assemblages ( Fig. 1d and Supplementary Figs. 2e and 4)-also began to shift in the late 1800s, followed by even more marked changes in the late-1960s (Breakpoint 1968). These trends are matched by accumulation rates of biogenic silica (BSi) (first breakpoint 1926; second 1978) (Fig. 1d), consistent with the overall increase in algal production inferred from the carbon proxies (C burial, C/N, and δ 13 C).
As the study lakes are shallow (mean depths = 2.7-4.9 m) and relatively clear (light attenuation, K a = 0.33-0.75 m −1 ) (Table 1), autotrophic production is dominated by benthic algae 31,32 . This conclusion is supported by related studies of nearby lakes, where benthic diatoms together with Nitella and mosses contributed substantially to whole-lake metabolism 33 . Whole-lake (pelagic + benthic) production estimates in these lakes ranged from 4 to 20 g C m −2 yr −1 (ref. 18 ).
Given that the organic sediment in our study lakes is predominantly autochthonous in origin (see ref. 34 for a discussion of the effect of benthic algae on isotopic signatures; see also Fig. 2a), the observed increases in C burial strongly imply greater aquatic production, which in turn requires increased nutrient input. Lakes in the Kuparuk Basin are generally Nlimited 35 , and nitrogen plays a central role in both terrestrial and aquatic productivity, as well as controlling community structure in this region 36 .
When interpreting sediment-core profiles, it is important, however, to recognise that their recent organic matter content is influenced to varying degrees by post-depositional diagenesis. Sediment C/N ratios initially reflect differential mineralisation rates, but as shown by Galman et al. 37 C/N profiles stabilise very quickly, within 10 years or so. Therefore, the emphasis here has been placed on longer-term changes. Bulk sedimentation rates have also increased (see ref. 38 ), which tends to counter diagenetic effects due to more rapid burial of organic matter. Furthermore, using a multi-proxy approach serves to emphasise the ecological signal as opposed to diagenetic noise. Similarly, the range of morphometries and water chemistry of the 5 lakes compared to the coherence of the geochemical response suggests that diagenesis is not the main process driving the sediment records. We therefore contend that the change in carbon burial in our study lakes is linked closely to changes in N dynamics in catchment soils, as well as inputs of reactive atmospheric N (Nr).
Nitrogen biogeochemistry. The change in sediment N-isotopic signatures (δ 15 N) (Fig. 1e) also began in the late-nineteenth century (Breakpoint = 1847-1878; Supplementary Fig. 3), which is temporally consistent with increased N loading to area lakes as reflected in the increased C accumulation rates, δ 13 C, and lake productivity indicators (Fig. 1c, d and Supplementary Fig. 2). As δ 15 N of bulk organic matter records complex changes in N fluxes from multiple sources as well as changes in N cycling, the direction and magnitude of the δ 15 N trends vary considerably among our five lakes ( Supplementary Fig. 2d). Specifically, two records show enrichment with 15 N (Perfect, Forgetful), one (Efficient) changes little, while two (Surprise, Relaxing) exhibit declining δ 15 N (Fig. 1e). Although the latter pattern is consistent with an influx of isotopically depleted Nr from the atmosphere 39 , alternate patterns in these and other regional lakes 40,41 demonstrate that N cycling is unlikely to be driven purely by atmospherically-derived N. The δ 15 N trends are unrelated in any significant way (Pearson's r = −0.48 to 0.57 with all p-values > 0.32) to physicochemical characteristics of the study lakes (Table 1), indicating importance of the differences in catchment N sources and in-lake N fractionation processes. Consistent with this interpretation, we note that baseline (~1800) δ 15 N values vary from sources characteristic of recently fixed N 2 (~0 ‰) to those in which terrestrial N is subject to intense recycling (~4 ‰) 42 .
That said, we note a striking pattern in the study lakes between the relative increase in sediment N content and the change in δ 15 N between the nineteenth century (1800-1880) and the last two decades represented in the cores . The lake cores with the largest increase in sediment N (ΔN) (Perfect, Forgetful) show a positive shift in δ 15 N, while the lakes with the smallest increase in ΔN (Relaxing, Surprise) exhibit a negative shift in  70 (green filled circles). Limnological response: c Mean lake-wide carbon burial in the five study lakes integrated by decade (grey fill); and z-scores for carbon: nitrogen ratios (open squares) and carbon isotopic composition (filled black circles) in master cores from the same five lakes. Fitted trend lines and 95% confidence intervals (shading) from a general additive model (GAM). d z-scores for biogenic silica (BSi) (open circles) and axis-1 scores from a PCA ordination of fossil diatom assemblages (green filled circles) for cores from two of the study lakes. Fitted lines and confidence intervals from GAMs. e z-scores for nitrogen isotopes in cores with increasing (Perfect, Forgetful; blue symbols) and decreasing (Surprise, Relaxing; red symbols) trends. GAM trend lines and confidence intervals as in Panels c and d. For the environmental response variables in the sediment-core records, two time-lines of major ecological change (vertical lines), the first beginning around 1880 (carbon burial and stable isotopes) and the second~1970 (diatoms and biogenic Si) are based on breakpoint analysis (see "Methods" and Supplementary information). δ 15 N over time (Fig. 2b). Although the exact mechanism(s) for this strong relationship (r 2 = 0.93) are uncertain, it suggests either differences in N-limitation (greater in lakes with the largest N increases) or in-lake N processing including N-fixation (greater in lakes with the smallest N increases), which is known to be important in Toolik area lakes 32 . Possible N sources are atmospheric Nr inputs as indicated elsewhere in the Arctic 39,43 (Fig. 1b), and perhaps, most importantly in moist tundra, lateral transfer of mineralised-N released from thawing soils 9,24 . Permafrost thaw may also increase export of terrestrial P to lakes, a process that is important to offset the P-limitation of primary production imposed by increased N influx.
Synthesis: coupled terrestrial-aquatic biogeochemistry. There is now unambiguous evidence for widespread environmental change in the Arctic 44 , especially from studies on the North Slope of Alaska. Greening (shrub expansion), thawing permafrost and increased thermokarst erosion have occurred in association with marked increases in regional air temperature since 1974 22,45 . The temporal disruption of C-N coupling within tundra landscapes, however, is not well constrained at decadal timescales 46 , such that longer-term evaluation requires a combination of ecosystem models 47 and palaeolimnological methods 48 .
Here, we have identified two main periods of change in regional biogeochemistry and C cycling (Fig. 1b, c): the first starting in the mid-to late-nineteenth century as indicated by rising C burial and declining δ 13 C values (breakpoint = 1877; Supplementary Fig. 3); and a second beginning in the midtwentieth century marked by accelerating C and BSi burial and changing algal composition (breakpoints for BSi 1926 and 1978, diatom turnover 1968). These temporal trends in sediment geochemistry are strikingly similar to those predicted by the results of ecosystem modelling 47 that showed terrestrial N mineralisation and loss in the Kuparuk Basin increasing after 1860 due to decreased soil moisture ( Fig. 1b; see "Methods") 47 . Typically, much of the N released from thawing soils is lost to the atmosphere through denitrification 24,49 , although some is also exported to streams and lakes where it can increase aquatic productivity 50,51 , as reflected in the sediment record.
Both terrestrial and aquatic ecosystems in the Kuparuk Basin are highly sensitive to environmental forcing 50 , and our palaeolimnological results suggest that biogeochemical cycling of C and N first began to change in the late-nineteenth century (Fig. 1c, e) well before anthropogenic warming of the atmospheric was substantial. Although a reversal of the circum-Arctic cooling associated with the Little Ice Age (LIA), which spanned the period approximately 1300 to 1850 AD, is well documented, most regional warming occurred after~1950 29 (Fig. 1a) due to a combination of human influence and a change in the Pacific Decadal Oscillation in the mid-1970s 52 . Further, the agreement between the two independent ecological hindcasts-modelled ecosystem N balance 47 and sediment geochemical recordssuggest that coupled aquatic-terrestrial C-N-dynamics in the Kuparuk Basin began to change in the late-nineteenth century driven by post-LIA changes in ecosystem and soil N processes (Fig. 1b) due to warming and altered soil moisture (see discussion in McKane et al. 47 ). The sediment-core results thus highlight both the protracted nature of ecosystem response to global-change forcing at high latitudes and the importance of a holistic approach to terrestrial-aquatic linkages during the Anthropocene, especially of altered N cycling 7,53 .

Methods
Study sites. Study sites are located on the North Slope of the Brooks Range in Alaska, USA, in the vicinity of the Toolik Lake LTER Field Station 38 . The local vegetation is dominated by tussock tundra, and there is continuous permafrost. The study lakes are small (<10 ha) and relatively shallow, with maximum depths of 5-13 m (Table 1). Lake waters are slightly alkaline (pH~8) and moderately dilute with conductivity ranging from 50 to 200 μS cm −1 . DOC ranges from 2.7 to 4.3 ppm. There are no distinct inflow streams, and catchment:lake area ratios are small (1.2 to 4.5). The present-day lake ice-free period in this area is~3.5 months, and monitoring data at nearby Toolik lake show no trend in either ice-off or ice-on over the last 20 years (https://toolik.alaska.edu/).
McKane et al. 47 suggested that N mineralisation in tundra soils of the North Slope of Alaska is strongly dependent on moisture content. They provide a detailed discussion of the role soil moisture and other factors influence the model results. Much of the nitrogen released from thawing soils is lost from the terrestrial ecosystems through denitrification, and terrestrial uptake is not able to utilise the expanding nutrient pool, such that some portion is transferred laterally (associated with increased hydrological flux) to streams and lakes where it complements atmospheric Nr deposition and in situ fixation thereby increasing aquatic productivity. The latter change is recorded in the lake-sediment record.
Field and laboratory methods. Details of core retrieval and core processing can be found in Fitzgerald et al. 38 . Briefly, three piston cores from each lake were fully dated using 210 Pb and analysed for organic content using standard loss-on-ignition methods 54 . 210 Pb chronologies were determined for the last 150-200 years using the constant rate of supply (c.r.s.) model 55 and extrapolated if needed to ca. 1800. 210 Pb dating accuracy declines substantially in sediments older than 150 years, but because of very slow sedimentation rates and near-constant supported 210 Pb, we were able to calculate dates for most cores back to c. 1800. The uncertainty for these dates is, of course, quite large, as shown in age-depth profiles in the Supplementary Fig. 1. Fig. 2 Contemporary aquatic biogeochemistry. a A comparison of carbon isotopic composition and carbon: nitrogen ratios in modern (core-top) sediment from the five study lakes with tussock peat (moss) and runoff water samples from the Toolik Lake catchment; and b a comparison of the change in sediment N concentrations with the change in nitrogen-isotopic composition between the nineteenth century (1800-1880) and the last 2 decades represented in the cores . The capital letters refer to the first letter of the five study lake names.
The dry mass sediment accumulation rate (DMAR, g cm −2 yr −1 ) was multiplied by the proportional biogenic silica (BSi) and organic matter (OM) content to derive BSi and OM accumulation rates. A correction of 0.469 54 was applied to the latter to derive the organic carbon accumulation rate (CAR; g C m −2 yr −1 ).
Samples of tussock peat (moss) and runoff water were collected from a 1-ha experimental sub-catchment of Toolik Lake in 2016 and analysed for C and N isotopes.
Whole-basin C fluxes. Carbon burial rates derived from single cores are not representative of whole-basin rates due to the problem of sediment focussing 56 . As multiple cores from each study lake were dated (see above and ref. 38 ), the cores could be integrated into a whole-basin rate and a focusing correction factor determined for each core. The focusing-corrected CAR were then integrated into decadal intervals and averaged for each lake. Reported modern organic CAR values exclude the most recent decade (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000), which may be biased high by incomplete mineralisation 57 .
Palaeolimnological proxies. A single, centrally-located core from each of the lakes was analysed for C and N isotopes, represented here in standard delta notation (δ 13 C and δ 15 N) relative to standard sources. Two additional proxies, biogenic silica (BSi) and fossil diatoms, were analysed on the central core for two of the lakes with good diatom preservation (Efficient and Perfect). Diatoms were poorly preserved in sediments of the other three lakes and so were not analysed. Biogenic silica content was determined using a timed-course alkaline extraction method 58,59 . Diatoms were identified and enumerated following standardised methods 60 , and community composition summarised as the first axis score in a PCA ordination. Stable isotopes were determined at the Institute of Environmental Change and Society, University of Regina, using a ThermoQuest, (F-MAT) Delta PLUS XL mass spectrometer coupled to a Carlo Erba elemental analyser and following Savage et al. 61 . δ 13 C values were not corrected for the Suess effect because the inlake CO 2 is mainly derived from recycled terrestrial carbon. The results of the Marine Biological Laboratory General Ecosystem Model (MBL-GEM) was used as an independent estimate of changing nitrogen dynamics in catchment soils and vegetation 47 .
Statistical analyses. Palaeolimnological data from each lake-sediment-core were transformed to z-scores for the purpose of inter-lake comparison, and the main changes in each proxy were undertaken on the mean z-scores of all lakes for a given variable ( Supplementary Fig. 3). The rationale for this approach was that it allows identification of departures from regional trends in the geochemical data through time across all or groups of lakes, while removing lake-specific variability. However, Supplementary Fig. 5a-c show breakpoints estimated on each core individually (untransformed) for those proxies that were z-transformed for the main breakpoint analysis (δ 13 C, δ 15 N, and BSi).
Major trends in palaeolimnological proxies were characterised by fitting generalised additive models (GAMs) 62 with a thin plate regression spline 63 . This method allows for non-linear modelling of trends with a formal assessment of uncertainty based on 95% confidence intervals 64 . To assess major areas of change in the geochemical and biological sediment records, breakpoint analyses were conducted to detect structural deviations from a simple linear regression model through time 65 . This method determines the optimal segmentation through an algorithm that computes the residual sum of squares for all possible break combinations given the minimum segment size of the time series (set at~5-7% of the observations, dependent on sample size). Ninety-five percent confidence intervals were calculated for the observations associated with each breakpoint date following Bai and Perron 65 . GAMs were modelled using the "mgcv" package and breakpoint analyses were conducted using the "strucchange" package in the R Statistical Software 62,66,67 .