Growing season leaf carbon:nitrogen dynamics in Arctic tundra vegetation from ground and Sentinel-2 observations reveal reallocation timing and upscaling potential

.


Introduction
Arctic tundra vegetation is undergoing substantial structural and functional change due to intensified climate warming in recent decades (Berner et al., 2020;Myers-Smith et al., 2020). Warming generally increases Arctic vegetation growth (Bjorkman et al., 2018a), thereby augmenting the relative importance of water and nutrient sources, primarily nitrogen (N) (Chapin et al., 2005;Shaver et al., 2006). In pristine Arctic tundra, plant available N is primarily a result of spatial reallocation, which occurs all the way from catchment scale (Skrzypek et al., 2015;Westergaard-Nielsen et al., 2020) to internal reallocation within individual plants (Evans, 2001). Processes within deciduous plants include both N uptake and transfer to branches and emerging leaves during spring as well as reallocation to stems and roots during autumn senescence. The latter process serves to store N during winter (Chapman et al., 2006), but it is not yet well quantified at the plant species level. Increased knowledge of N availability and temporal patterns in N cycling and allocation strategies in dominant tundra plant species is therefore important for understanding growth feedbacks in Arctic ecosystems resulting from a warmer climate.
The ratio between carbon (C) and N (C:N) in plant leaves is a proven tool to assess and understand plant N dynamics and resource optimization strategies (Mack et al., 2004;De Deyn et al., 2008;Semenchuk et al., 2015). However, there are a lack of studies focusing on species-specific leaf C:N dynamics with sufficient temporal resolution and extent to pinpoint important periods of changes. Equally important, we lack insight and tools to estimate the net effects of plant N-dynamics over larger tundra areas (Beamish et al., 2020;Myers-Smith et al., 2020), which is otherwise vital for our ability to project impacts of future ecosystem feedbacks to climate.
Warming-induced changes in tundra plant community composition can affect carbon-climate feedbacks in several ways. Observations and experimental open-top chamber (OTC) manipulations show that increased plant growth leads to greater above-and below-ground C storage (Elmendorf et al., 2012a;Elmendorf et al., 2012b;Sistla et al., 2013;Zamin et al., 2014). However, plants are not just responsible for C uptake, they also determine soil microbial decomposition rates, and thereby soil C storage, through species-specific differences in functional traits, such as leaf area and nutrient content (Bjorkman et al., 2018a). Moreover, differences in canopy height and leaf size, which again is linked to nutrient availability, affect ecosystem energy balance and soil thaw depth (Chapin et al., 2005;Blok et al., 2010). This has implications for microbial activity through enhanced regional warming or localized shading-induced cooling, respectively. As foliar plant litter is the most prominent input of C to the soil (Aber and Melillo, 2001), changes in the chemical composition of litter fall regulates nutrient cycling rates and soil organic matter storage (Christiansen et al., 2018a(Christiansen et al., , 2018bParker et al., 2018). Low leaf C:N ratios, and high cellulose:lignin ratios, have long been known to catalyse decomposition (Melillo et al., 1989). Consequently, interactions and feedbacks between litter fall, soil N mineralization, and plant C dynamics in a climate change perspective have been increasingly studied at local Arctic sites for decades (e.g. van Heerwaarden et al., 2003;Jonasson et al., 2006;Semenchuk et al., 2015;(Blok et al., 2016) Christiansen et al., 2018aChristiansen et al., , 2018b. Despite recent efforts in synthesising circumpolar-wide ecological responses to environmental change (Bjorkman et al., 2018a;Prevéy et al., 2019), it remains unclear how different tundra plant species will respond to ongoing warming in terms of build-up of plant biomass C and N. While single site studies are of great importance for improving our process understanding, it is challenging to upscale the findings from local sites to landscape-and regional scales for quantifications of the aggregated climate feedback effects. This is particularly important when considering recent evidence showing complex plant responses to climate warming at local scales (Myers-Smith et al., 2015;Phoenix and Bjerke, 2016;Myers-Smith et al., 2020). Here, one of the grand challenges is our ability to assess detailed nutrient contents and status of tundra canopies across space and time with regards to C:N ratios.
Remote sensing approaches facilitate effective monitoring of biomewide tundra change. Vegetation indices from satellite remote sensing (such as NDVI) associated with plant cover have been used to quantify Arctic tundra aboveground biomass (Jia et al., 2006;Raynolds et al., 2012;Westergaard-Nielsen et al., 2015;Räsänen et al., 2019), and for descriptive phenology during the Arctic growing season (Zeng et al., 2011;Karlsen et al., 2014;Park et al., 2016;Karami et al., 2017). Using satellite data, Park et al. (2016) reported an advancement in the onset of the Arctic growing season of 2.4 days per decade since 1982, oftentimes resulting from earlier snowmelt dates (Zeng et al., 2011;Westergaard-Nielsen et al., 2017;Pedersen et al., 2018). While optical satellite-based remote sensing can provide important data to indicate ecosystem changes and are scalable to global coverage, satellite imagery alone provides limited information on how climate warming continues to shape remote Arctic plant communities. Thus, challenges remain when using remote sensing tools to quantify plant biophysical or biochemical properties, such as plant height or leaf nutrient status (Beamish et al., 2020;Myers-Smith et al., 2020). In addition, satellite imagery alone does not shed light on functional phenology, e.g. what impact does climate change have on seasonal variations in canopy nutrient status and -reallocation, specific parameters causing plant stress, or absolute photosynthesis rates. For that kind of functional detail, we need on-theground observations. This is an ongoing issue because in-situ validation datasets with high temporal and/or spatial sampling frequency remain rare for many Arctic regions, due to the challenging logistics associated with field campaigns in this part of the world. Nevertheless, recent plotscale observations suggest diverging vegetation responses to climate change across the Arctic, e.g. with increased growth and primary productivity in some areas and decreased in others (Phoenix and Bjerke, 2016;Vowles and Björk, 2019;Myers-Smith et al., 2020). These responses are highly difficult to capture with e.g. NDVI (Myers-Smith et al., 2020), emphasizing an important need for linking large-scale remote sensing methods to detailed vegetation dynamics on the ground.
Efforts have been made to link different spatial scales with focus on measuring vegetation qualities linked to their functioning, e.g. canopy N, C, and lignin. These efforts typically require either access to hyperspectral data (Martin and Aber, 1997;Smith et al., 2002;Martin et al., 2008;Hollinger et al., 2009;Knyazikhin et al., 2013) which constrains the spatiotemporal resolution and scalability, or include advanced modelling of the canopy radiative transfer system (Boegh et al., 2013;Jay et al., 2017) which rely on high-quality input datasets. However, literature reporting the performance of these methods to estimate dynamics of e.g. canopy N, C, or lignin in Arctic ecosystems is to our knowledge non-existent. This has led to an urgent need and requests for tools to better estimate functional plant changes from space (Myers-Smith et al., 2020), based on accessible satellite data with high spatiotemporal coverage, and validation with ground observations. In an Arctic context this lack of knowledge also stands in contrast to the relevance for herbivores such as reindeer (Rangifer tarandus) which rely on both plant quantity and feeding quality, and to the current trend of increased farming potential at high northern latitudes following climate warming. This pinpoints an urgent need to explore and optimize grazing and farming practices in the Arctic (King et al., 2018a(King et al., , 2018bHannah et al., 2020).
Here, to the best of our knowledge (Beamish et al., 2020), we present the first growing season-long quantification of temporal fluctuations in Arctic tundra leaf C:N in combination with a remotely sensed scaling approach. We do this by pairing in-situ plant species-specific observations, obtained approximately weekly from June through October in Arctic Greenland, with imagery using a satellite-derived spectral index. This approach allows us to assess the empirical spatiotemporal relationship between in situ leaf C:N contents from four of the most common tundra shrub species and the integrated plant community-based spectral signal achievable from the 10-20 m spatial resolution S-2 satellites. In short, we (1) present and validate a novel methodology that is scalable from local plot to regional scales for dwarf shrub heath vegetation that is widespread across the circumpolar Arctic and (2) provide a unique quantification of plant species-specific dynamics of N, including estimates of when and how much of total above-ground leaf N mass is reallocated from leaves to stems and roots during autumn. This novel method advances our remote sensing capabilities to include assessments of temporal functional changes in tundra vegetation concerning leaf N content, timing of internal plant N-redistribution, and litter quality during senescence. Consequently, this study has implications for our understanding of climate-driven changes in Arctic tundra vegetation, regarding plant communities and C fixation as well as ecosystems with grazing, including future ambitions of human-introduced domestic grazers in the Arctic.

Methods
Our study was carried out from June through October 2016 in the Østerlien valley (69.25359 o N, 53.51301 o W), located on the southern part of Disko Island in West Greenland (Fig. 1). The in-situ study area covers approximately 0.1 km 2 on a gentle south-facing slope going from 60 m.asl. to a few m.asl towards the Disko Bay shore. The Østerlien valley lies just east of the research station Arctic Station, and the area A. Westergaard-Nielsen et al. surrounding Østerlien has been studied systematically since 1996. Arctic Station's environmental monitoring programme became part of the Greenland Ecosystem Monitoring programme in 2018 (https://g-e-m. dk/).

Site climatology
The study area lies in a discontinuous permafrost zone (Obu et al., 2019) with mean annual summer (June-August) and winter (December-February) air temperatures of − 3 ± 1.8 • C, 6.8 ± 1.3 • C, and − 12.5 ± 3.9 • C, respectively (1991-2017 means ± SD; Zhang et al., 2019). From 1991 to 2017, the mean annual air temperature increased significantly by 1.6 • C per decade, mainly due to a pronounced winter warming of 3.2 • C per decade. A marked warming trend occurred in 1991-2005, while data for 2006-2017 show greater year-to-year variation with no significant trends in mean annual, summer, or winter air temperatures which is in line with general temperature trends in the region , see Zhang et al. (2019) for further details. Mean annual precipitation is estimated at ~400 mm (Hansen et al., 2006). The soil is well-drained and poorly developed from weathered basaltic bedrock parent material.

Site vegetation
Mesic dwarf shrub heath tundra dominates the Østerlien valley. The vegetation consists mainly of deciduous shrubs Betula nana, Salix glauca, and Vaccinium uliginosum, and the evergreen shrub Empetrum nigrum. Some graminoids, suchs as Carex sp. and Poa sp., and various mosses and lichens are present but to a much lesser extent than the dominant shrubs. This dwarf shrub dominated vegetation type belongs to the Circumpolar Arctic Vegetation Map (CAVM) subunits P1, P2, and S1, making the study area representative of up to 25% of the ice-free circumpolar low Arctic (CAVM subzones D and E according to Walker et al., 2005), including substantial parts of the ice-free Greenland .

Experimental setup and field sampling
In early June 2016, we established nine 1 × 1 m plots following a distributed sampling scheme with a minimum distance of 65 m between plots. We established each plot so that it provided representation of the vegetation in a surrounding area of 20 × 20 m to allow subsequent comparison with S-2 satellite data. Based on spectral (red, green, blue) UAV data at 3 cm spatial resolution, we found the 20 × 20 m area around each plot to have a near gaussian distribution, and that the mean pixel value from the plots was within one standard deviation of all pixel values in the 20 × 20 m area. From this we concluded that selected plots represented as homogeneous an area as possible. The nine plots were dominated by different combinations of vascular plants Betula nana, Salix glauca, Vaccinium uliginosum, and Empetrum nigrum. We visited the plots 13 times during the 2016 growing season, with a median of eight days between visits from June 14th through October 15th. During each visit, we photographed the plots at nadir approximately 1.5 m above the surface (G1X camera, Sony, Tokyo, Japan). Each image covered an area of ~3 m 2 . In addition, we randomly sampled 10 leaves of each of the four dominant plant species mentioned above. Leaves were collected at random from plant individuals growing within 1 m of each plot. Note that no leaves were collected inside plots. This approach allowed us to sample leaves continuously throughout the growing season, minimizing sampling stress on specific individuals, while obtaining unbiased photographs from within each plot.
We used photographs from 20th of July to estimate the percentage cover of each of the four dominant plant species by averaging a visual estimate of species cover from a 6 × 6 grid drawn across the photograph (supplemental fig. S1). Photographs from the remaining dates served for visual inspection of disturbances and phenology. The combined cover of the four shrub species was above 88% in all plots so we standardized cover of these four species to 100% cover, which also aided in compensating for leaf N in other, non-sampled, plant species. The remaining surface cover in the plots consisted of a combination exposed bedrock, mosses, lichens, and graminoids. The site coordinates were logged 5 times (GPSmap 62, Garmin Ltd., Kansas, USA), and the Euclidian averages were selected as the representative coordinates to decrease spatial uncertainty.
In addition to leaf sampling, site-level soil pore water was sampled 12 times at the Met Station ( Fig. 1) during the study period for subsequent quantification of dissolved organic nitrogen (DON), NO 3 − and NH 4 + . The soil nutrient sampling procedure and analysis followed Rasmussen et al. (2020).

Leaf carbon and nitrogen analyses
The sampled leaves were air-dried and stored at ambient conditions approximately one year before undergoing lab analyses at the Center for Permafrost (Copenhagen, Denmark). In this study, we consider the petiole as being a functional part of the stem (Pérez-Harguindeguy et al., 2013), and therefore to have no marked influence on the remotely sensed reflectance compared to the leaf blades. Consequently, we carefully removed petioles from the leaves in the lab. Subsequently, we homogenized (Retsch MM200 ball mill, Retsch Gmbh, Haan, Germany) five leaves from each plant species at each of the 13 sampling dates, whereafter 3-4 mg subsamples were analyzed for C and N content in an elemental analyzer (Flash 2000, Thermo Scientific, Bremen, Germany) coupled in continuous flow mode to an isotope mass spectrometer (Thermo Delta V Advantage, Thermo Scientific, Bremen, Germany). We used pure gases of CO 2 and N 2 calibrated against certified materials of 15 N-(NH 4 ) 2 SO 4 and 13 C-sucrose, as certified reference materials.
Species-specific leaf C:N ratios and N content were subsequently derived as an average of all the analyzed leaves per species across the nine plots per sampling time. We report all leaf data as community-weighted means per plot, using the standardized percent cover of the four dominating plant species (Table 1), or we report data as species-specific means across all plots.
We used plant species-specific leaf mass per area (LMA) data, obtained from the Tundra Trait Team database (Bjorkman et al., 2018b) and calculated as specific leaf area (SLA) − 1 , in order calculate canopy biomass N for each of our plots and sampling times. We used circumpolar-wide medians of LMA (using means instead of medians We used database LMA data in combination with our own plot-scale leaf N content data and vegetation cover to calculate communityweighted peak-summer canopy N mass and autumn senescence N mass as follows: Note that our estimation of Canopy N does not include N in aboveground stems. We did not estimate Canopy N for dates prior to July as leaves are still developing in June and our calculation used database LMA values that were collected primarily in July and August, i.e. when leaves are near to fully developed. Our estimation of Canopy N assumes that LMA remains constant during autumn senescence, which likely underestimates resorption efficiency and overestimates remaining N mass in near-senesced leaves (van Heerwaarden et al., 2003). Cover fractions of individual plant species varied between plots (Table 1), and by deriving a standardized weighing-factor for each species in each plot, we amplify the absolute importance of the present species in subsequent estimates of canopy N. The cover standardization builds on the assumption that the plots are fully covered by vegetation and that other plant species, which were not accounted for in our analyses (e.g., forbs or graminoids), have a similar N content to the shrubs investigated by us. While these relative cover estimates could add bias to our in-situ leaf C:N, the four sampled shrub species always covered >88% of the plots, and only Plot 9 had minor barren patches (approximately 3% of the area). Moreover, other plot-scale plant community cover analyses show that forbs and graminoids each cover ≤10% of the area in Østerlien (Christiansen, unpublished data), with forb and graminoid leaf %N being broadly similar to our three sampled deciduous shrub species (data from the Tundra Trait team database show that those species range between 1 and 4% N by leaf dry mass). We therefore argue that the alternative to standardizing the four measured species to a 100% cover (i.e. to assume no C or N in the plant species that are unaccounted for in our analyses but still constitute a minor part of the surface cover) is suboptimal and would lead to underestimation of canopy N. We acknowledge that our plant cover estimates based on top-down plot orthophotos likely underestimate total leaf abundance due to our limited ability to account for layered canopies (vertical variation and overlapping canopies/leaves). In Østerlien, all four sampled shrub species usually grow within multi-layered canopies, and our canopy N mass calculations do not account for N in sub-canopy leaves.

Satellite imagery
Satellite data were acquired from the sun-synchronous S-2A satellite, featuring the Multi Spectral Instrument with 13 spectral bands, and a spatial resolution between 10 and 60 m (Drusch et al., 2012). Due to converging satellite paths at higher northern latitudes, the satellite revisit time can be as low as 3 days in our study region, The image data were downloaded through the Copernicus Open Access Hub, and all imagery was radiometrically and atmospherically corrected using the Sen2Cor module with a maritime atmosphere setting (Louis et al., 2016). The satellite scenes in supporting table ST1 were found to match our criteria of maximum 50% cloud cover and proximity in time to the insitu sampling dates (+/− 2.3 days on average). We then extracted cloud-free satellite pixels, covering our field plots for subsequent comparisons. Band 6 (red-edge, 740 nm) and 11 (shortwave-infrared (SWIR), 1610 nm) were used to compute a normalized reflectance index (NRI) at 20 × 20 m spatial resolution: where ρ is bottom of atmosphere reflectance at 740 and 1610 nm, respectively. The index is based on a correlogram between the 13 satellite bands and in-situ C:N values, from which bands 6 and 11 (at 20 m spatial resolution) showed among the highest correlation coefficients (Pearson's rho). In addition, our approach is motivated by the nitrogen SWIR absorption features around 1510 and 1690 nm (Herrmann et al., 2010;Dechant et al., 2017), and the well-documented ability of the rededge region to describe vegetation type and status (Horler et al., 1983;Gitelson et al., 2006;Meneghin et al., 2018). In addition, the red-edge is associated with chlorophyll, and thus indirectly with N content (Clevers and Gitelson, 2013). Although band 11 is not centered at either absorption wavelength, it overlaps with 1690 nm. Note that the SWIR bands, as well as the near-infrared bands are also affected by moisture levels on the Earth's surface. This feature is e.g. exploited using a Normalized Difference Wetness Index (NDWI) as proxy for soil moisture (Jørgensen et al., 2015) or plant canopy moisture (Gu et al., 2008). Other normalized difference indices combining SWIR and red-edge bands have been used to estimate vegetation water stress and canopy water content, e.g. Fensholt and Sandholt (2003) who used the MODIS band 2 and 6, or Zhang et al. (2018) who used S-2 band 8a and 11 to derive normalized difference infrared indices (NDII). Therefore, we computed NDWI based on S-2 band 12 and 8 (Gu et al., 2008), three NDII indices based on S-2 bands 8a and 11, bands 7 and 11, and bands 8 and 11 as well as the commonly used NDVI based on the standard band 4 and 8 (Erinjery et al., 2018), to test the potential for isolation of unwanted noise from changes in moisture regime.
The S-2 satellites can suffer from inconsistent spatial accuracy at high northern latitudes due to inconsistencies in the elevation model used, which can vary with satellite orbit (Kääb et al., 2016). In our study area, we found shifts to be <10 m between orbits (satellite scenes were visually compared to a sub-meter accuracy orthophoto), possibly resulting from our study area being relatively flat, and consequently our in-situ plots were still located in the same pixels regardless of orbit. We were able to confirm this theory also by using the constant visible location of the nearby Arctic Station.

Data analyses
We used R v. 4.0.0 (RCore, 2016) to process data and calculate all canopy N estimates. Unless stated otherwise, we used least squares linear regression and t-tests for the correlation analyses (using Matlab 2018a, Mathworks, Natick USA).
When observing data points following an increasing or decreasing trend over time, we can estimate when in time a certain ratio of change in the data from their min to max has occurred. To quantify this timing in our satellite and in-situ datasets, we normalized the data to [0-1], interpolated linearly between data points, and extracted the dates when an increase or decrease (represented by the interpolation model) had reached 0.4, 0.5, and 0.6 of its change, respectively (using Matlab   Table 1 Standardized cover fraction distribution of the sampled shrub species in nine replicate plots in Østerlien valley, Disko Island, West Greenland. Prior to standardization, the total relative cover of the four shrub species was >88% for all plots. 2018a, Mathworks, Natick USA).

Plot scale plant cover
We estimated cover fractions of the dominant shrub species based on camera photographs from each plot and its nearest surroundings as defined by the photograph's coverage. S. glauca was present in all plots, while V. uliginosum and B. nana were the second most frequent plant species (Table 1).
From visual inspection of weekly photographs, the relative canopy leaf area per species did not change noticeably from July to late September. Leaf coloring shifted from green to carotenoid-dominated between the 3rd and 10th of August for plot 2, 3, and 6, and between the 10th and 24th of August for the remaining plots.

Local growing season climate, soil water chemistry, and satellite NDVI
The rapid drop in soil moisture in early July reflects the shift from early-season saturated soils, due to snowmelt and subsequent lateral runoff, to a mesic soil moisture regime in summer ( Fig. 2A). Plot-scale NDVI, derived from the S-2A data, followed an expected pattern with a rapid green up in June, following snowmelt, followed by a senescing phase starting in mid-August. Visual shifts in vegetation color, based on plot photographs, appeared long before the first freezing temperatures around the 20th of September. Leaf color change followed the period of decreasing temperatures starting from late August ( Fig. 2A and supplemental S2), which was also reflected as the period with significant increase in NRI 1610 for all plots. This temperature drop coincides with 3 precipitation events on the 28th of August, 8th of September, and 30th September respectively, seen as local increases in soil moisture at 10 cm depth ( Fig. 2A).
The summer of 2016 was largely an average summer in terms of soil moisture as well as air and soil temperatures according to Zhang et al. (2019), although spring snow cover was relatively low due to a warm early spring. This partly resulted in a drier summer than normal.
Soil water samples taken during the study period were analyzed for DON, NO 3 − and NH 4 + . Throughout summer, the total mass of the organic N fraction were factors higher than inorganic N chemically bonded in NO 3 − and NH 4 + in the soil water (Fig. 2B). However, all three components followed a similar seasonal pattern with relatively high availability in early spring, followed by a significant decline (below our detection limit for both solute concentrations and total mass of NO 3 − and NH 4 + ) during peak growing season. DON, and to some extent NH 4 + , increased markedly in late summer until freeze-in in early October, presumably due to leaching of litter fall. Solute concentrations are reported in supplemental fig. S3.

Seasonal dynamics in leaf C:N ratios
The analyzed leaf samples showed an expected difference in C:N through time between the deciduous and evergreen species. Specifically, E. nigrum had higher leaf C:N in the early and peak growing season than the three studied deciduous species, mainly due to lower N levels ( Fig. 3A-C or supplemental fig. S4). In contrast, deciduous shrubs B. nana, S. glauca, and V. uliginosum showed a moderate increase in C:N starting already from the first sample on June 14th, and with a more distinct increase in the first half of August until the end of sampling on October 15th (Fig. 3A). For B. nana, we saw an unexpected decrease in C: N on the last sample date, however the difference falls within the standard error margin, and lab tests showed that the drop in C:N was caused by variation between the random leaves picked for analysis at this late phenological stage.
At the plot level, the shrub species cover varied according to Table 1. In general, the cover of S. glauca and V. uliginosum was most prominent across plots. The temporal change of leaf C:N at plot level is summarized in Fig. 3A. Looking at the N fraction only, S. glauca leaves contained about 2 / 3 of the leaf N in August compared to July, and by early September, S. glauca leaves contained 1 / 3 of peak summer %N. B. nana leaves reallocated leaf N slightly earlier than S. glauca as 2 / 3 and 1 / 3 of the highest measured leaf N remained in leaves by 12th of July and 22nd of August, respectively. For V. uliginosum, the corresponding dates were 29th of July and 5th of September, respectively. For the evergreen E. nigrum, leaf N dynamics followed a different pattern relative to the deciduous shrubs. E. nigrum leaf N content remained stable throughout summer, except for a slight decrease starting from late August (Fig. 3B).

Seasonal dynamics in canopy N mass
We calculated seasonal dynamics in canopy N mass for each shrub species (not shown) and for each plot based on in-situ leaf %N (Fig. 3B). The three deciduous shrub species showed a consequent decline in canopy N throughout the growing season, while canopy N for the evergreen E. nigrum was more stable over time, which is also reflected in measured leaf %N. When looking at community-weighted means at the plot-scale, plot canopy N varied between plots due to differences in species cover fraction. Plot 2 and 3 showed a distinct drop in mid-to late August, roughly halving their N. Plot 6 had an even earlier significant drop from early August, while the remaining plots displayed later drops, or a quite linear decline through the season (Fig. 4).

Results from analyses of satellite indices
A regression analysis between satellite-indices from plot 5 (nearest to the soil moisture point measurement) showed that the near-surface soil moisture at 10 cm ( Fig. 2A) and 30 cm depth was not correlated with NDWI as defined in section 2.6, nor with NRI 1610 , NDVI, or the NDIIs. All tested satellite indices showed a weak correlation with soil moisture at 20 cm depth, and show thereby the potential importance of soil characteristics for the temporal variation in indices. However, we note that the correlation arise from data distributing as one period with stable satellite indices yet markedly decreasing soil moisture (July), and a subsequent period with marked change in satellite indices and stable soil moisture (August-September), suggesting a need for a dedicated experimental study setup to investigate the dynamics (supplemental fig.  S5). Other soil environmental characteristics may be similarly important but are beyond the scope of this study.
We found that satellite-derived NDVI from the nine plots was significantly negatively correlated with the in-situ C:N ratios through the studied season, using one regression model per plot (Table 2). When evaluating the NRI 1610 , again as one regression model per plot, we found a significantly better correlation (Fisher's Z transformation, p = 0.037) with in-situ C:N. A regression analysis comparing in-situ C:N and absolute values of NRI 1610 (supplemental fig. S4) for all plots and sample times in one single model showed a significant correlation, however with visible uncertainty around the mean values of C:N and NRI 1610 . This uncertainty greatly diminished when comparing C:N with a normalized NRI 1610 dataset, yielding an 81% level of explanation (Fig. 5), and suggesting that the slope of linear fit lines varied between plots. Substituting NRI 1610 with absolute NDVI showed a significant relationship, as did using a min-max normalized NDVI signal, so the explanatory level of the normalized NRI 1610 was significantly better than NDVI (Fisher's Z transformation, p = 0.04). Two apparent outliers on Fig. 5B represent plot 3 and 6 showing a local earlier rise in C:N compared to the other plots (Fig. 5B). NRI 1610 showed the highest degree of explanation of in-situ C:N compared to the three other computed NDIIs, but the difference was not statistically significant (NDII results in supplemental table ST2).
When studying the temporal evolution of NRI 1610 during the growing season we observed a small increase on average for the nine vegetation plots from the beginning of the sample period until mid-August, followed by a substantial increase from mid-August, and a leveling-off in mid-September. In absolute numbers, Plot 6 and 9 had a significantly   higher initial NRI 1610 values in the early season than the other plots, while also reaching the highest index values in late autumn (supplemental fig. S4A -C). A min-max normalization per vegetation plot allowed for a direct comparison of the range of shift over time between plots, showing lower variation between plots (Fig. 5B).
Since the general patterns between in-situ leaf C:N and NRI 1610 were in good agreement, we extracted transitional dates from the temporal trends based on min-max normalized data of in-situ C:N, NRI 1610 as well as NDVI. When using a 50% threshold, i.e., the date at which half of the seasonal change had occurred, we saw a statistically significant linear correlation where the timing of a 50% change in NRI 1610 could explain 82% of the timing for 50% change in C:N. The correlation was significant at 40 and 60% seasonal change as well (Table 3), although with greater residuals (see supplemental fig. S6 for scatter plots). The NDII based on band 7 and 11 showed an equal R 2 to NRI 1610 (R 2 = 0.82) but a marginally higher mean square error of 16.06 (see supplemental table ST3 for all NDII results). There were no significant relationships with C: N transition dates when using satellite-based min-max normalized NDVI.
Using the correlation between NRI 1610 and leaf C:N (i.e., one regression model for all plots) we derived three maps, covering the studied heath area (Fig. 6).

Discussion
Plant uptake of available N and interspecific strategies for autumn Nreallocation are critical N-cycling pathways in terrestrial Arctic ecosystems. Differences in the amount and timing of internal N-reallocation will influence how distinct shrub species cope with current and future climatic changes, with implications for herbivore forage quality and availability as well as ecosystem C and nutrient cycling processes. Yet, methods to assess large-scale temporal and spatial dynamics in canopy N contents are lacking (Beamish et al., 2020). Here, we demonstrate growing season and autumn variation in leaf C:N ratios at a relatively high level of temporal and spatial resolution, ensuring capturing of N reallocation patterns at the plant species level as well as integrated plant community effects for mixed dwarf shrub heath tundra. These findings are based on a development in remote sensing-based methods that allows us to take the first steps towards an upscaling of seasonal canopy C: N dynamics for large-scale assessments of functional plant phenology and leaf N dynamics in Arctic tundra.

Scaling leaf carbon to nitrogen patterns from plot to satellite imagery
By including S-2 satellite data, we target the knowledge gap between in-situ measurements of N-dynamics at the small-scale leaf process-level and estimates of impacts at larger landscape spatial scales. However, linking across spatial scales raises the question of representativeness. The sampled heath type is the second-most widespread vegetation class on Disko Island, covering ~600 km 2 . On a pan-Arctic scale, this vegetation class corresponds to prostrate-shrub tundra (as defined by Walker et al., 2005), covering up to 25% of the icefree land area across the Arctic (Walker et al., 2005). Being able to capture plant functional dynamics in this type of dwarf shrub heath from space would therefore be substantial for our understanding of climate feedbacks related to changes in vegetation composition (Loranty et al., 2014), spatial dynamics in carbon fixation and N-cycling in vegetation (Myers-Smith et al., 2020), and the perspective of increased extensive livestock grazing at high northern latitudes (King et al., 2018a(King et al., , 2018b. Using a surface classification map by Karami et al. (2018) and our estimation of total canopy N, and assuming that autumn reallocation starts from August 10th (see Fig. 5A), we find a local reallocation of ~0.8 g N m − 2 through the senescing phase before litter fall (8.7 mg N m − 2 day − 1 ). This resorption corresponds to roughly 50% of peak growing season canopy N, similar to leaf N resorption efficiencies of deciduous tundra plant species (Chapin III et al., 1980;Chapin III and Moilanen, 1991) Scaling to dwarf shrub heath on the Disko Island, this corresponds to a biological mechanism likely transferring ~480 t N from leaves to stems and roots, and the remains to the soil surface to be recirculated through microbial breakdown of leaf litter.
The NRI 1610 showed great potential for explaining the general temporal variation of C:N across plots, and it was statistically significantly better than NDVI. NDVI did not correlate with the derived C:N transition  Table 3 Results of transition date timing correlation analyses between satellite indices and in-situ leaf C:N (mse = mean square error). dates at which a certain level of C:N change (presumably due to N reallocation) had occurred. In contrast, the NRI 1610 functioned as a highly promising proxy for identification of these transition dates. Of the other three tested NDIIs, NRI 1610 performed only marginally better than when using band 7 and 11, showing that the use of a red-edge band is important to capture C:N dynamics, likely as a combined effect of vegetation type and the link between chlorophylls and N. Our findings therefore support the use of spectral indices other than NDVI. The use of SWIR wavelengths in NRI 1610 questions to which extent the observed spectral changes result from changes in functional plant phenology or changes in surface moisture. Other studies have confirmed and utilized the SWIR sensitivity to estimate surface moisture at high northern latitudes (Bartsch et al., 2009;Jørgensen et al., 2015), however, most often as a proxy for conditions in the top surface 0-5 cm. (Sadeghi et al., 2017). In our case, there was no correlation between soil moisture at 10 cm and SWIR, nor did we observe a visual covariance over time (supplemental fig. S5), which could otherwise be used to disentangle the spectral signal into a moisture and plant-related component. By normalizing the NRI 1610 signal, we are confident that we have mitigated any major potential biases between plots resulting from systematic differences in moisture regime. Plant species vary in reflectance within the red-edge portion of reflected light, however, within-plant temporal dynamics such as changing leaf chlorophyll/carotenoid relationships also affect reflectance (Gitelson et al., 2006;Meneghin et al., 2018). Although the red-edge portion is barely visible to the human eye, the visible earlier shift in autumn coloring in plot 2, 3, and 6, relative to other plots, could be indicative of a relatively earlier N reallocation, even if there was visible variation in leaf color within plots. Fig. 5A shows an earlier marked increase in C:N for plot 3 and 6, suggesting a visible proxy for this phenomena. This is also evident on the NRI 1610 signal for those plots (Figs. 5B) which shows a similar earlier index increase, although it does not fully mimic the range of C:N shift as seen from the two outliers in Fig. 5. We assign this to within-plot variation in leaf senescence stages, illustrating the challenge of representative leaf sampling during autumn senescence. Nevertheless, due to the close association between N and chloroplasts (Roberts et al., 2012), part of these mechanisms are arguably captured by the red-edge spectrum.
Finding and calibrating remotely sensed spatial proxies for the leaf N-dynamics is an important step towards bridging gaps between in-situ observations and ways to estimate the actual dynamics and implications over a larger spatio-temporal domain. The correlation between leaf C:N changes and NRI 1610 does not prove a direct link between N-content in leaves and their spectral properties. Indeed, the plants will transport and reallocate other important nutritional elements alongside N, e.g., phosphorus (P) and potassium (Chapin III et al., 1980), which could influence the observed spectral changes. However, what we show here for the first time is indeed a directly scalable satellite-based index which correlates with C:N reallocation mechanisms in a widespread tundra vegetation class where N is a critical limiting factor for growth (Weintraub and Schimel, 2005;Semenchuk et al., 2015;Pedersen et al., 2020), allowing for novel spatial assessments of plant reallocation strategies. Similarly, our approach describes the timing of a key functional phenological event in its autumn senescing phase, which is oftentimes more complex to capture using optical remote sensing than e.g. spring greenup (Zhang et al., 2018).

Temporal dynamics in species-specific leaf carbon and nitrogen content
As expected, we observed distinct growing season differences in leaf C:N dynamics between deciduous and evergreen shrub species. The deciduous shrubs all contained relatively high concentrations of leaf N at the first sampling time in mid-June, with leaf N declining throughout the growing season. This pattern indicates that the movement of N from seasonal storage in roots and stems and into new leaves during the first couple of weeks following snowmelt (i.e. at bud-break) had finished by our first sampling (Chapin III et al., 1980). During this time in June, soil soluble N availability peaks to a level likely controlled by climate, before rapidly declining to near-zero for the remainder of the growing season (Fig. 3C), indicating N limitation (Semenchuk et al., 2015;Rasmussen et al., 2020). Therefore, as leaves matured and increased in size during June and July (supplemental fig. S3), leaf C:N ratios gradually increased for all three deciduous shrubs due to competition with soil microbes for a decreasing amounts of soil N and a corresponding dilution of N with recent photosynthesized C compounds (Fig. 3A). Leaf C:N ratios began to increase at much greater pace in August, due to reallocation of N from leaves to stems and roots for winter storage, with the onset and rate of N reallocation being slightly offset between the deciduous species. B. nana lost half of its max leaf %N a week earlier (~25th of July) than V. uliginosum and two weeks before S. glauca, respectively. This is partly due to a higher initial leaf %N in B. nana than for S. glauca and V. uliginosum. The time when one-fifth of max B. nana leaf N remained was again a week earlier than V. uliginosum and two weeks before S. glauca, respectively. Overall, this may suggest an earlier reallocation timing for B. nana and, to a lesser extent, V. uliginosum compared to S. glauca, although with some spatial variability between plots. Species exhibiting a slightly more opportunistic strategy N-use in late autumn are interesting in the light of increased relative climate warming in autumn, as has been served in Greenland  and in the Arctic in general (Walsh, 2014;Bintanja and Andry, 2017), although further research is needed to explore such plantspecies/climate feedbacks. In contrast to deciduous shrubs, tundra evergreens typically have less biomass belowground and they tend to keep nutrients stored within their canopy (Chapin III et al., 1980), as exhibited here by relatively stable E. nigrum leaf N, and only subtle C:N change, throughout the measurement period. Fig. 6 illustrate a first step towards spatiotemporal mapping of C:N, and cover a relatively homogenous area of tundra heath. Generally, the whole site exhibits an increase in C:N from peak season (July) to senescence (October), in line with individual satellite cells extracted from the nine plots. The highest spatial variation seems to occur during the early phase of senescence (September), in line with the between-plot variation observed in Fig. 5B. From a visual inspection it seems the 20 m spatial resolution could be sufficient to describe variation at a hectarescale, however, there are indications of spatial autocorrelation, which we interpret as mainly resulting from autocorrelation in plant species composition. Using a single equation for all plots to model C:N does not capture the true variation. The calibrated model (linearly calibrated to in-situ C:N) yields much improved absolute results, underlining the critical importance of ground validation.
The patterns of N forms in the analyzed soil water samples showed similar dynamics to previously reported seasonal changes, with relatively higher availability in early spring, followed by a rapid decline likely resulting from increased microbial activity and plant photosynthesis (Roberts et al., 2009), and a near-zero plateau during the peak growing season. The rapid drop in soil water N during spring and summer indicate a critical N limitation (Semenchuk et al., 2015;Rasmussen et al., 2020) underlining the relevance of using vegetation C:N dynamics as an indicator of N location in the general ecosystem. Lastly, we measured a brief increase in NH 4 + and distinct increase in DON in the autumn. Such increases have previously been attributed to autumn freeze-thaw cycles damaging roots and thus limited plant N uptake (Christiansen et al., 2012). However, in our case it is more likely due to increased N input from litterfall (Semenchuk et al., 2015;Blok et al., 2018) which is further supported by the occurrence of rain events in September-October increasing the mobilization of ions from litter to soil.

Implications for natural and domesticated grazers
Seasonal changes in leaf nutrient concentrations have significant consequences for herbivore food preference. Arctic tundra ecosystems are nutrient-limited and plant species differ in tolerance levels and strategies to cope with nutrient availability. Even minor changes in soil nutrient pools may therefore alter competitive fitness of certain plant species. This can lead to plant species-specific winners and losers in a climate-change context, which may have major implications for both natural and domesticated grazers. For example, increased shrub growth may alter forage quantity and quality for caribou (Rangifer tarandus) (Zamin et al., 2017). In summer, Arctic ungulates selectively forage deciduous shrubs (Boertje, 1984;Manseau et al., 1996), such as Betula spp. and Salix spp., due to their common abundance and relatively high nutrient content compared to other shrubs (Mårell et al., 2006;Skarin et al., 2008). Changes in canopy N content are particularly important as greater leaf N provides greater nutritional value for herbivores. Consequently, warming-induced shrub dominance is likely to impact wild and domesticated R. tarandus populations through changes in forage quantity and quality (Zamin et al., 2017). Here, phenological changes are also very likely to affect herbivores in the Arctic.
Leaf nutrient contents are generally high when leaves emerge in spring, whereas nutrient contents decrease continuously until senescence in autumn, as seen here in this study and in others (Chapin III and Moilanen, 1991). For migrating herbivores, such as geese, early snowmelt leads to earlier onset of graminoid growth, and may result in decreased leaf N content when geese arrive (Doiron et al., 2015). As it is critically important for herbivores to maximize nutritional gain during the short Arctic growing season, trophic mismatches therefore pose a severe threat for herbivore population fitness during this decade (Couturier et al., 2009). This is also an issue for animal husbandry, mainly sheep in South and West Greenland, which relies on the quality of natural vegetation during summer when animals are grazing open lands. Here, optimal forage quality and the timing of moving animals indoors is a critical parameter, which can be optimized by improved knowledge of canopy nutrient content. This is particularly relevant when assessing year to year variations before investing in infrastructure such as sprinkler systems to allow for longer seasons, increased number of animals per area, or even mowing (King et al., 2018a(King et al., , 2018bHannah et al., 2020). The methodology and results presented here show how such scaling for larger regions can be made based on on-going and free remote sensing products. We consider results presented here as the first step towards considerably improved quantifications of spatiotemporal variations in tundra canopy N content (or C:N ratios). Several critical knowledge gaps exist such as the performance in wetlands, and a decoupling with soil moisture regime, however it is a promising advancement for our efforts to manage both natural and domesticated herbivore populations now and in the future.

Implications for future sampling procedures, satellite-based studies and modelling of carbon dynamics
The marked changes in leaf C:N ratios reported here, at the plot community and individual plant species levels during the entire growing season, highlight that reported leaf C:N ratios should be related to a specific time of the growing season. Single snapshot measurements are not representative for the plants in question, not even relatively, as slightly divergent plant-specific reallocation strategies can occur. It also confirms the importance of picking leaves for litter degradation studies at the latest possible stage to avoid a bias towards too high litter Nlevels, which in turn may result in overestimated breakdown rates and therefore underestimated C sequestration rates.
We propose that future studies should focus on experimental treatment plots with nutrient additions, e.g. N or P fertilizer addition studies, and possibly include greater spectral and spatial resolution image data in order to explore direct links between spectral properties and element cycling within tundra plants further. Such experimental plots need to be of sufficient size to accommodate the spectral resolution of available satellite imagery, and while this might be challenging there are certainly sites with plot sizes large enough to accommodate the needs. Moreover, it is arguably important to know the approximate cover fraction of species or similar diversity measures such as remote sensing-derived leaf area index estimations. It is also worth noting that our findings only relate to a specific vegetation class, and we see an urgent need for including contrasting climates, year to year variations, and different vegetation types to allow for calibrated scaling across even wider regions. Lastly, this study focuses on C:N ratios and future work should include and address the importance and scaling of other critical plant leaf characteristics, including protein content and lignin which is essential for plant quality as a forage source.

Conclusions
This study is an important first step towards bridging the gap between mechanisms of N-cycling at specific plant species and plot community levels in Arctic tundra with the integrated landscape results observed spectrally from space. The eligibility of the study is thereby as A. Westergaard-Nielsen et al. a tool for quantifying spatiotemporal heterogeneity in dwarf shrub heath canopy C:N ratios. Similarly important, it demonstrates the importance of season-long sampling strategies for quantifying speciesspecific leaf N dynamics, both to improve our process-level understanding of tundra N-cycling and to guide future sampling strategies and modelling efforts linking leaf N contents to soil C sequestration rates. We recommend further controlled experiments to establish direct links between spectral properties compatible with easily accessible satellite data and element cycling, and we encourage further calibration efforts to extrapolate the usage of this new spectral methodology for other vegetation classes and moisture regimes.

Author contribution
Bo Elberling (BE), Casper Tai Christiansen (CTC), and Andreas Westergaard-Nielsen (AWN) designed the study and layout of field work. CTC conducted in-situ leaf sampling, plot photography, and canopy N modelling. AWN was responsible for analyzing leaf and satellite/UAV data as well as climate data. All authors contributed to writing the manuscript.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.