Common mechanisms explain nitrogen-dependent growth of Arctic shrubs over three decades despite heterogeneous trends and declines in soil nitrogen availability

 Heterogeneity has been observed in the responses of Arctic shrubs to climate variability over recent decades, which may reflect landscape-scale variability in belowground resources. At a northern fringe of tall shrub expansion (Yuribei, Yamal Peninsula, Russia), we sought to determine the mechanisms relating nitrogen (N) limitation to shrub growth over decadal time.  We analysed the ratio of 15 N to 14 N isotopes in wood rings of ten Salix lanata L. individuals (399 measurements) to reconstruct annual point-based bioavailable N between 1980 – 2013. We applied a model-fitting/model-selection approach with a suite of competing ecological models to assess the most-likely mechanisms that explain each shrub’s individual time-series.  Shrub δ 15 N time-series indicated declining (7 shrubs), increasing (2 shrubs) and no trend (1 shrub) in N availability. The most appropriate model for all shrubs included N-dependent growth of linear rather than saturating form. Inclusion of plant-soil feedbacks better explained ring width and δ 15 N for eight of ten individuals.  Although N trajectories were individualistic, common mechanisms of varying strength confirmed the N-dependency of shrub growth. The linear mechanism may reflect intense scavenging of scarce N; the importance of plant-soil feedbacks suggests that shrubs subvert the microbial bottleneck by actively controlling their environment. and (ii) overall inter-series correlation Both statistics calculated


Accepted Article
This article is protected by copyright. All rights reserved The ecological mechanisms driving shrub-nitrogen relations may operate on varying temporal resolutions, including longer timescales than current observations. In order to predict future dynamics of shrubs in the Arctic, we need to understand the potential for N availability to constrain growth over decadal time, which is currently poorly represented in land system models (Bouskill et al. 2014). Space-for-time substitution analyses (in which environmental gradients of soil N and shrub growth are measured to indicate how processes may operate through time) overestimate the rates and order of processes (Elmendorf et al. 2015), owing to each location's varying history. Without observations over decades, it is unclear if short-term processes are masking longer-term N limits to growth (Luo et al. 2004).
Research has demonstrated the use of stable N isotopes in tree rings as an integrated indicator of pointbased soil N availability (McLauchlan et al. 2007;Gerhart & McLauchlan 2014). Here, we adapt this approach to Arctic shrub wood rings, which can be precisely dated using dendrochronological techniques (Forbes et al. 2010;Myers-Smith, Hallinger et al. 2015); their slow growth however poses challenges to the N isotope method through low annual wood production. By combining ring width and N isotope data from shrub rings, it becomes possible to reconstruct individual-based time series of local soil N availability and plant secondary growth. This method (McLauchlan et al. 2018) is extremely useful as it allows us to infer past plant-N relations beyond the spatial and temporal constraints of long-term experimental plots.
Here, we sought to demonstrate a method to assess the nature and strength of N limitation to shrub growth covering the period shrubification in recent decades. To achieve this, we measured annual time-series of stable N isotopes from shrub wood rings covering a 35-year period that matches the remote sensing record (Tucker et al. 2005), and is of sufficient length to capture processes of plant-N relations beyond short term physiological responses to variability in soil N. We measured individual time-series from ten individuals of the same species -Salix lanata L. -from a site on central Yamal Peninsula (Error! Reference source not found.) totalling 399 individual N isotope measurements, to establish inter-shrub variability in absolute values of the N isotope proxy, taking into account the number of individuals (often five or less) used for oxygen and carbon wood ring isotope proxies (McCarroll & Loader 2004). We used the time-series to test four key hypotheses: H1. Soil N availability has been increasing over the previous three decades H2. Secondary growth rates of Salix lanata L. shrubs have been N-limited by the rate of root foraging for N H3. Plant-soil feedbacks are important to explain the rate of Salix lanata L. secondary growth H4. Tall shrubs are nearing geometrical limits to stem diameter

Accepted Article
This article is protected by copyright. All rights reserved

Site
Our study site was Yuribei, Yamal Peninsula, Russia (68.91 °N,70.23 °E). Between 1980 and 2019 the Marre Sale weather station (69.72 °N, 66.82 °E -161km northwest, 24m elevation) indicates that the mean July temperature was 12.4℃ (95% between 8.5℃ and 16.2℃), and mean January temperature was -23.4℃ (95% between -29.9℃ and -15.1℃); the wettest month was August (48.22mm precipitation) and driest March (12.95mm precipitation). Temperatures rise above zero for only four months of the year (June-September), although since 1980 summer temperatures have warmed 1℃ to 2℃ (Menne et al. 2012). The soil is a Histic Cryosol characterised by poor drainage (Mikhajlov 2017). Yuribei is classified as low-shrub tundra (Subzone D) in the Circumpolar Arctic Vegetation Atlas (Walker et al. 2009). The plant community consists of deciduous shrubs and graminoids. Salix lanata L is the dominant shrub species in lowland areas: in recent decades this species has increased substantially in height to upwards of two metres (Macias-Fauria et al. 2012).
During a field campaign in July 2013, 52 shrub discs were collected from the riparian lowland at Yuribei (Error! Reference source not found.). Shrub stems were sampled to obtain basal discs at the closest level to the ground. Only larger individuals were sampled for the purposes of ring-width chronology development. Effort was taken to sample independent genets; when possible, neighboring female and male plants were tested for potential variation in stem allocation. We measured ring width increments on wood discs to a precision of 0.01mm. Standard dendrochronological techniques were employed to age and crossdate each individual shrub disc, precisely identifying missing rings (Fritts 1976). Thus, each ring measurement is dated to a calendar year (see Methods S1). Relative increment time-series were calculated as the current year's ring increment divided by the cumulative radius. We also collected data in Central Yamal to calibrate allometric equations (see Section 3.3.4). Shrub size metrics were measured on 200 Salix lanata individuals: these included stem length, stem basal circumference, and stem diameter.

N Time Series Development
δ 15 N in bulk wood has previously been validated as a proxy of point-based soil nitrogen availability (McLauchlan et al. 2017;Gerhart & McLauchlan 2014;McLauchlan et al. 2007). From our 52 shrub samples, we selected ten at random from which we developed time-series of δ 15 N using the dated annual wood increments (Error! Reference source not found.). Of these, three (shrubs YUSL03, YUSL05, and YUSL09) grew on slightly elevated gentle slopes (i.e., within 500 metres of an elevation greater than 15

Accepted Article
This article is protected by copyright. All rights reserved metres above sea level (a.s.l.)), whereas seven were sampled from the lowest elevations in the landscape (i.e., flat plains below 10 metres a.s.l.), growing on thickets close to the Yuribei river main and adjacent channels. We conducted a t-test to establish if the mean of all δ 15 N measurements differed between the two elevations (within 500 metres of elevation >15 metres a.s.l. vs further than 500 metres distant from >15 metres a.s.l.). A major methodological consideration for wood δ 15 N chronologies is temporal resolution: as the quantity of bulk wood required is high (~11mg) in comparison to δ 18 O and δ 13 C wood isotope measurements, pooling of annual increments into bins from increment cores is usually required to obtain the necessary wood (Gerhart & McLauchlan 2014). However, we sought to develop annual time-series by taking advantage of the higher wood quantity available from wood discs versus increment cores.
The shrub discs did not have symmetrical growth forms: for each shrub disc, a 15mm-wide wood block was sectioned following the longest radial from bark to pith. The height of wood blocks varied between 8 and 15mm. To increase visibility of the annual ring boundaries, the surface and reverse of the wood block were wet with de-ionised water for 10 minutes, then smoothed using a GSL1-microtome (Gärtner et al. 2014). A thin section was mounted as a cutting guide using standard wood anatomical methods (Gärtner & Schweingruber 2013) (see Method S2). Each wood block was separated into annual increments using a Stanley blade on a granite cutting surface under a Leica DSM 300 macroscope. Cuts were made in the direction of the grain and verified against ring boundaries on both surfaces of the wood block to minimise sampling error between adjacent annual increments. Bulk wood from each individual increment was homogenised by: (a) finely chopping the sample using a Stanley blade, and (b) milling the sample into a wood powder using a Spex MixerMill 8000M. Milled samples are more homogenised thus reducing standard deviation of isotope measurements (Riechelmann et al. 2014). Each binned sample was milled in a steel capsule using two ball-bearings for 60 minutes. From each resultant powder, 11.00±0.05mg was weighed into a 5x8mm tin capsule. Replicates were weighed for one in three samples to characterise sample variability. The δ 15 N and %N content of each sample was analysed in a random order with a Thermo Fisher Delta V+ mass spectrometer coupled to a Carlo Erba NC 2500 Elemental Analyzer at the Central Appalachian Stable Isotope Facility, University of Maryland Centre for Environmental Science.
We identified significant trends in δ 15 N time-series by fitting linear models based on Theil-Sen single median using the mblm R package (Komsta 2013) under R 3.6.3 (R Core Team 2020). The Theil-Sen statistic is especially suited to short and noisy time-series (Hoaglin et al. 2000). We conducted an intercomparison of time-series to ascertain if there were common trends: we calculated (i) mean between-shrub correlation ( ), and (ii) overall inter-series correlation (in which a mean chronology is built serially through stepwise addition of shrub individuals). Both statistics were calculated using methods from Bunn (2008).

Accepted Article
This article is protected by copyright. All rights reserved

Modelling of Plant -Nitrogen Interactions
To assess if shrub growth was N-limited, the role of feedback mechanisms, and if shrubs are nearing geometrical limits to stem diameter, we used a Model-Fitting and Model-Selection (MFMS) approach (Bonsall & Hastings 2004) to infer the underlying mechanisms operating on shrub growth and nitrogen availability. MFMS approaches compete alternative ecological models to determine the most appropriate model structure to explain the observed time-series, while also estimating the strength of the mechanisms involved. We constructed a set of mechanistic top-down models (Poorter et al. 2013) that relate shrub secondary growth to a single limiting resource (Error! Reference source not found.). Within the modelling approach, we addressed specific considerations that enabled the use of dendroecological data: (a) allometric relations between biomass and basal diameter / ring width (Section 3.3.4); and (b) use of δ 15 N as a proxy for N availability (Section 3.3.5).
We adapted Tilman's resource-ratio models (Tilman 1990;Tilman 1988;Huisman 1994;Jabot & Pottier 2012) to capture the main processes in plant-nitrogen dependence through continuous-time ordinary differential equations (ODEs). We assume that a shrub's ability to grow is a function of its photosynthetic capacity, minus respiration costs, and loss of existing material to the environment (e.g. via litter). Given that N is an essential macronutrient (Tilman 1990), we used a nitrogen-economy approach (Garnier 1991) where nitrogen availability may limit the production of new biomass as a single limiting resource. Soil nitrogen availability was assumed a function of background nitrogen replenishment in the soil, minus the nutrient requirement of any plant material grown, minus environmental losses relative to the amount of nitrogen present. Here, background N replenishment included inputs from biotic processes (i.e. decomposition, N mineralisation, N fixation) and abiotic processes (i.e. N becoming available from thawing permafrost, N deposition), which are outlined in Error! Reference source not found..

Change in Nitrogen = Background Replenishment -Plant Uptake -Environmental Loss [+ Plant Inputs]
The model was defined by a base model, with three substitutable functions representing (a) N limitation, (b) plant-soil feedback, and (c) geometric constraints. In the base model, shrub biomass, B, and soil nitrogen dynamics, N, was represented by a dual-equation ODE system:

Accepted Article
This article is protected by copyright. All rights reserved Here, is a resource-constrained net growth function (which accounts for photosynthesis and ( ) respiration); is a geometric effect on the growth rate; represents an intrinsic growth rate of plant ( ) biomass (i.e., a growth constant representing plant performance); and represents the rate of environmental loss of biomass from other factors (e.g. litter). is a linear background replenishment rate of soil nitrogen; is a plant-soil feedback mechanism that depends on current biomass; and is a ( ) density-dependent loss term for nitrogen from soils. We included two additional key assumptions to extend this approach to wood ring data. First, allometric relations determine the proportionality of basal diameter to total biomass. Second, the basal diameter of the main stem can never decrease once grown, whereas shrub mass may decline. The individual components, the product of which form the model hypotheses, are outlined in Table 1.

Mode of Resource Limitation f(N)
To identify if shrub secondary growth has been limited by N uptake, we included three forms of Nlimitation: growth where N is not limiting; a linear resource-dependence; and a resource-dependence based on root foraging. Within these models, we conceptualised nitrogen limitation as occurring via two mechanisms: (a) limits on photosynthetic rate; and (b) limits to the shrub's ability to scavenge for N. First, low N availability can be linked to reductions in leaf health; it has previously been deduced that, below a saturation point, tissue N concentration relates to photosynthetic productivity with a Monod saturating function (Van Rees 1994), where K is a half-saturating constant. Second, photosynthetic demand for N may be greater than the efficiency of roots to scavenge for and uptake available N. Root uptake rate is monotonically increasing, and usually also represented by a Monod or Michaelis-Menten function (Jabot & Pottier 2012), requiring a similar constant K that represents the halfsaturation N required for growth at half of the maximum rate, and a scalar a that represents the N-use efficiency: McNickle and Brown (2014)

Accepted Article
This article is protected by copyright. All rights reserved Here, represents N-use efficiency, and h represents the 'handling' or processing time for each unit of N (McNickle & Brown 2014). Assuming that plant root systems act as a colony of individual 'foragers' for N, the individual forager nitrogen uptake rate is multiplied by the units of root biomass, . Given that both photosynthesis and uptake may be rate-limited by N, and both can share a Holling disc functional form, the operation of both processes is here abstracted to a single function: Where gross biomass production (i.e, without losses) = and shrub N uptake = . Here represents the combined effect of nitrogen uptake and conversion into plant biomass. The parameter r represents the intrinsic growth rate (i.e,, the efficiency of converting N into biomass), a represents N-use efficiency (i.e., the efficiency of root uptake), and h represents the rate at which the combined processes of uptake and incorporation occur. This implies that given the similarity in the form of the modelled functions between both processes, a saturating model in our study would not distinguish between photosynthetic rate limitation or rate-limited root uptake, but would clearly indicate the efficiency of each process individually and the levels of N availability able to saturate one (or both) of the two ecological processes.

Presence of Plant-Soil Feedbacks
To identify if plant-soil feedbacks were important determinants of growth and soil N, we included an optional feedback loop. Salix lanata have ectomycorrhizal associations in Yamal (Akhmetzhanova et al.

2012)
, which may access organic N forms to overcome in part the microbial bottleneck (Chapman et al. 2006). We represented plant-soil feedback as an addition to the N stock, where the addition is proportional to the current biomass loss rate. Thus, the strength of the feedback depends on current shrub biomass. We assume that the relation is linear, following Jeffers et al. (2012): Here, represents a conversion factor of plant biomass into soil N, and represents environmental loss as a fraction of total plant biomass .

Accepted Article
This article is protected by copyright. All rights reserved

Presence of Geometric constraints
To identify if tall shrubs are nearing geometrical limits to stem diameter, we included an asymptotic growth function where there is a limit to shrub biomass. In dendrochronological studies of Arctic shrubs, ( ) age-related / ontogenetic growth trends are often absent (e.g. Macias-Fauria et al. (2012)). Assuming that the maximum size of woody plants is constrained by mechanical limits, this indicates that either: (a) such shrubs are far from their asymptotic stem diameter and growing slowly, thus appearing to grow linearly over the observed timeframe; or (b) shrubs preferentiate investment in new young stems versus old stems, thus never approaching asymptotic size. We represented mechanical constraints using a Richards' form, as this non-linear plant growth model has a mechanistic basis (Paine et al. 2012): Here, B is current plant biomass, and K represents asymptotic plant mass owing to mechanical or other constraints.

Allometry
To gain a mechanistic understanding of the processes driving shrub ring production, it is necessary to apply an allometric model; ring width is a one-dimensional size variable of a plant stem, compared to twodimensional cross-sectional area, and three-dimensional volume (Burkhart & Tomé 2012). As the conceptual models used here require information on total biomass, we apply an allometric model to quantify the size-specific three-dimensional stem volume laid down for each millimetre of one-dimensional wood increment. We connected ring width to stem mass using a shrub allometric model that mimics fractal woody structure within a multi-stemming plant individual (Götmark et al. 2016

Accepted Article
This article is protected by copyright. All rights reserved

Nitrogen Availability
We applied an explicit transform of δ 15 N to standardised N availability using an existing empirical relationship. δ 15 N is an index relative to air (0‰)  . Working with the assumption that fractionation between δ 15 N in leaves and wood is insignificant, we transformed our δ 15 N data with the above global model. We also conducted a sensitivity analysis to understand the robustness of our modelling approach to the assumed globally-calibrated relationship between unitless soil N availability and δ 15 N, versus an alternative approach where the relationship is calibrated employing a curve obtained from subset of data from ectomycorrhizal species in temperate south-east Alaska and which, despite obvious major differences to our study system, was deemed to be the most appropriate localised subset of available data. Inter-site variability may arise from variability in local soil conditions, mycorrhizal types, and stages of soil development (Hobbie et al. 1999). However, we consider that the above global model is most appropriate for our analysis, as (a) the greater geographical, species, and mycorrhizal variability dampens noise arising from site specificity, and (b) a subset of data specific to the climate and Quaternary history of the Western Siberian Arctic is not available. A full sensitivity analysis for the global calibration versus the most appropriate local calibration -which uses data from Alaska -is illustrated in Supplementary Notes S1.

Accepted Article
This article is protected by copyright. All rights reserved

Model Fitting and Model Selection
We determined the most appropriate combination of N-limitation, plant-soil feedback, and geometrical limit to explain the ring width and δ 15 N time-series data for each shrub individual by applying a Model-Fitting and Model-Selection (MFMS) approach. First, we conducted maximum likelihood estimation (MLE) for each model in Table 1. MLE identifies the best set of model parameters given the data via an optimisation routine that maximises the goodness-of-fit. To assess goodness-of-fit, a negative log likelihood to be minimised was defined following a bivariate normal distribution: Here, n is the time-series length; is the observed cumulative stem radius; the observed δ 15 N at time t; the cumulative stem radius derived from biomass (estimated in the dynamical ODE model) then  Table S1). Starting bound ranges were specified based on exploration of the model properties; the ranges specified could generate time-series that reflected the magnitude of the observed time-series. The global optimisation routine used minimises the importance of starting bounds by allowing wide exploration of parameter space. 68% and 95% confidence intervals were estimated for each model parameter using a likelihood interval approach (Morgan 2000).
Model selection was conducted using constrained Akaike Information Criterion (AICc), through which the likelihood is corrected for the number of observations in short time series and penalised for the complexity of mathematical models. We used the AICc to calculate Akaike weights for each model and for each shrub individual (Wagenmakers & Farrell 2004). Akaike weights translate the AICc into a relative measure of

Accepted Article
This article is protected by copyright. All rights reserved evidence in support of each model given the data. In addition, we compared the goodness-of-fit of the shrub-nitrogen models between individuals by generating one-step-ahead predictions for stem radius and δ 15 N at each models' MLE, from which we calculated the root mean square error (RMSE) values of both cumulative stem radius and δ 15 N.

Dendroecological Data
Replicate δ 15 N measurements for one in every three years to determine variability within wood ring δ 15 N indicated an overall standard deviation of δ 15 N values of 0.303‰ (3 d.p.). The variability was high within and between shrub discs (Error! Reference source not found.). Where a shrub exhibited missing growth rings, we could not measure δ 15 N; this resulted overall in 11 missing measurements within the ten timeseries. The associated wood N concentration time-series are shown in Figure S2.

Trends and Variability in δ 15 N Time-Series
Nine of the ten δ 15 N time-series exhibited significant linear temporal trends (1980 -2013). A declining δ 15 N trend between 1980-2013 was found for seven time-series, while increasing trends were found in two time-series. For the seven declining series, the median decline was -0.036‰ per year -representing a median depletion in δ 15 N of 1.188‰ -between 1980 and 2013. When all ten δ 15 N were included, the overall median depletion of δ 15 N was 0.759‰. Theil-Sen summary statistics for all ten individual timeseries are shown in Table 2.

Accepted Article
This article is protected by copyright. All rights reserved We found spatial heterogeneity in the variability in δ 15 N. Mean inter-series correlation was 0.207 (3 d.p.) and mean between-shrub correlation was 0.055 (3 decimal places (d.p.)), suggesting that N availability was heterogeneous between the locations of shrub individuals. It was therefore not appropriate to construct a mean site-level chronology of δ 15 N. Despite the low overall inter-series correlation, some pairs of series showed higher correlations; was 0.46 between five of ten individuals. Landscape position was a factor in the spatial differences of absolute δ 15 N values, as higher mean δ 15 N values observed in the shrub individuals growing on the slope (YUSL03, YUSL05, and YUSL09) versus lower values in those collected from the lowland area (t-test: p = 0.0001).

Calibration of Allometric Relations
We determined the relationship between shrub height and stem diameter by applying a non-linear leastsquares fitting procedure to the shrub stem length -basal diameter dataset, using the Niklas and Spatz (2004)  with 99% significance (p<0.001). When this empirical stem length -basal diameter function was incorporated into the shrub radius -biomass allometric function, biomass increased exponentially with stem radius (Error! Reference source not found.).

Mechanistic Model Fitting and Selection
We conducted model-fitting and model-selection to identify the best model formulation to explain the observed δ 15 N and ring-width data for N-limitation, plant-soil feedbacks, and geometric effects. The Akaike weights for all models are shown in Table 3 (for likelihood values see Table S3).

Form and Strength of N Limitation
Akaike weights -which represent the relative support for each hypothesis -indicated that N-dependent functional forms (N2 and N3) better represented the data than N-independent growth (N1) for all ten shrub individuals assessed. The linear form of N-limitation (N2) was more appropriate than the root foraging (saturating -N3) form for nine of ten individuals. However, the strength of this mechanism was variable between individuals. Given the most likely estimated values of a, which represents the efficiency of N uptake, and their 95% confidence intervals, the variability was greater between individuals, with little overlap of the 95% confidence intervals (Error! Reference source not found.). For one shrub individual (YUSL39), the root foraging N-dependent growth function (N3) was found to better explain the observed data; the mean observed and predicted δ 15 N values were lowest for this individual.

Accepted Article
This article is protected by copyright. All rights reserved

Plant-Soil Feedback
The inclusion of a plant-soil feedback mechanism (F2) improved the ability of models to explain the observed data for eight of ten shrub individuals according to the Akaike weights. Models without a plantsoil feedback (F1) were the most appropriate for YUSL30 and YUSL39. For the eight shrubs that exhibited important plant-soil feedbacks, the feedback between plants as estimated by the litter conversion factor , which represents the N content derived from each gram of biomass, and the biomass loss rate , showed variable strength between the individual shrubs (Error! Reference source not found.). Plant-soil feedbacks (i.e., the inferred rate of process y(B)) contributed substantially to total soil N input (i.e., background N replenishment plus plant-soil feedback) over time compared to background N replenishment from other mechanisms' inputs (as estimated within the parameter ).

Geometrical Constraints to Stem Diameter
The most appropriate model included geometrical constraints (i.e., asymptotic limits) to maximum biomass (R2) for six of the ten shrubs, indicating that for these individuals stem growth is approaching geometrical limits. For four individuals -YUSL3, YUSL21, YUSL29, and YUSL34 -a linear functional form (R1) was more appropriate. Of the six shrubs for which a limit to maximum size was important, the median asymptotic basal stem diameter was 11.66cm, ranging between 7.36cm -38.15cm (all 2 d.p.). The equivalent median biomass was 22.03kg.

Goodness-of-Fit of Models to Data
The best-fitting models for all shrub individuals were able to replicate the observed trajectories of δ 15 N.
Model performance, however, was relatively poor for YUSL14, the only individual for which there was no significant trend in δ 15 N over time (Error! Reference source not found.). As an additional indicator of goodness-of-fit, we calculated one-step-ahead predictions for the best-fitting model for each shrub individual and calculated the root mean square error (RMSE) for δ 15 N and stem predictions. The RMSEs for the most appropriate models are shown in Table 4 (and for all models in Tables S4 and S5).
The RMSE of one-step-ahead predictions indicate predictive accuracy of the model in the original data units and reflecting the differences in variability between each time series. The median RMSE of δ 15 N for all ten individuals was 0.53‰. Similarly, the median accuracy for predicting stem radius was 0.72mm (see

Accepted Article
This article is protected by copyright. All rights reserved Figure S3 for predicted time-series). In addition, examination of the squared residuals between observed and expected stem radius indicated that model performance was worse during years in which the mean ring width chronology indicated relatively slow growth (i.e. 1988, 1989, and 1997).

Trends in Soil Nitrogen Availability
Our wood δ 15 N time-series indicate that soil N availability declined from 1980 to 2013 at multiple locations at Yuribei, contrary to our hypothesis. However, the presence of positive trends in two sampled individuals demonstrates heterogeneity of trends in N availability at the landscape scale. The presence of both declining and increasing δ 15 N trajectories suggests that the observed N availability has derived from localscale factors (e.g. local plant structure and microtopography) rather than regional atmospheric patterns, from which we would expect a common response. Given the rapidity of the warming that has occurred on the Yamal Peninsula over the recent decades, it is surprising that we found declining N signals over decadal

Accepted Article
This article is protected by copyright. All rights reserved 'terrestrial oligotrophication' (Craine et al. 2018). It is plausible that increasing the lengthening growing season in the Arctic may be responsible for similar increases in N demand.
Differences in absolute values of δ 15 N between individual time-series appear to confirm high spatial heterogeneity of soil N availability within the landscape (or 'patchiness' -localised differences in fertility arising from varying importance of specific N cycling processes) . Two lines of evidence -greater mean

Role of N in shrub growth
A key result is that we identified consistency of N-limitation and its mechanism for these Salix lanata individuals, although processes varied in strength. This suggests that common non-linear processes have governed relations between Salix lanata shrubs and variability in soil N within the period of recent climatic change. The importance of N-limitation at Yuribei aligns with evidence from tundra fertilisation experiments, which indicate increases in above-ground biomass production following addition of NPK fertilisers (Bouskill et al. 2014). That a linear rather than saturating N-limitation function better represented the data suggests that all bioavailable N is utilised during an individual growing season, with no evidence that N-limitation is tapering off. Thus, any limiting rate on N uptake and N incorporation into biomass (i.e. root scavenging and photosynthetic rate) occurs at a finer temporal resolution than our data. Indeed, as the which could indicate reliance on organic N sources.
The model-inferred rates and magnitudes of soil N transfers suggest that N inputs and losses are spatially heterogeneous at a fine scale. N-limitation to Arctic shrubs may be moderated by three major mechanisms:

Accepted Article
This article is protected by copyright. All rights reserved decomposition rate, competition with N-immobilising microbes, and access to resources from thawing permafrost (Wild et al. 2018). The identified dominance of plant-soil feedbacks as part of total soil N inputs over time (Error! Reference source not found.-Right) may be explained by increased access to organic deposits within thawing permafrost or biological N 2 fixation. As was relatively high on slopes, it suggests that such processes that govern N replenishment vary with topography. In addition, modelinferred rates of soil N loss -which may include gaseous losses, leaching, and microbial N immobilisation-varied substantially between the individuals studied, which may reflect fine spatial variability in loss processes. Indeed, evidence from active layer samples in the Siberian Arctic suggests that plant N limitation did not result from soil organic matter ( high, but unrelated to the incorporation efficiency of N into biomass (r): this suggests that shrubs here may respond to near-future increases in N availability by increasing biomass production without subsequent increases in litter quality, while also being unconstrained by N-uptake rates.

Shrubs as drivers of N cycling
The importance of shrub-soil feedback mechanism for all of the shrubs tested indicates that Salix lanata shrubs drive changes in their abiotic environment over decadal timescales. Given that the median estimated value of N-uptake efficiency (a) increased by an order of magnitude when shrub-soil feedbacks were incorporated into models, our findings corroborate with other evidence: Buckeridge et al. (2010) found that -in low Arctic tundra -soil N cycling rates were up to three times faster on tall-birch tundra than dwarf shrub tundra due to litter feedbacks. Our modelling results suggest that shrubs exert significant influence on local N-cycling rates. The identified variability in the strength of the plant-soil feedback between individual plants, which is reflected in the parameter values of B and ⍺ (Error! Reference source not found.), may

Accepted Article
This article is protected by copyright. All rights reserved indicate both the environmental circumstances of the individuals but also phenotypic plasticity in terms of litter quantity and litter quality.
Two lines of evidence -the importance of litter feedbacks, coupled with the linear form of N-limitationsuggest that the Salix lanata individuals subvert the microbial bottleneck by actively controlling their environment (Chapman et al. 2006). We suggest that Salix lanata are able to overcome low N conditions by accessing organic N forms via their ECM associates, as this would enable fast turnover of shrub litter.
Organic N sources can be significant for N nutrition in tundra environments: Hobbie and Hobbie (2006) applied a mass balance approach that quantified plant-mycorrhizal fluxes based on 15 N fractionation at Toolik Lake, Alaska, identifying that 61 -86% of N in plant biomass was obtained from mycorrhizal sources.

Variable Importance of Geometrical Limits
Given that the inclusion of asymptotes to growth only improved model fits for six of ten shrub individuals, it appears that some shrubs are nearing geometrical limits to stem diameter and that this is not caused by Nlimitation. For the shrubs where an asymptote was not included our results suggest that the role of geometrical limits was of lesser importance than N availability in controlling maximum size, or that these plants were far from their maximum size. In comparison, common plant growth models assume that the asymptote reflects a combination of mechanical and nutrient limits to growth (Paine et al. 2012). That the importance of an asymptote was unrelated to the shrub stem size during the modelling timeframe suggests that unseen ecological or behavioural (e.g. stem abandonment) mechanisms are driving stem thickening.
High variability in modelled asymptotic sizes between individuals suggests that mechanical limits (beyond which the shrub form cannot support itself) are not responsible; shading and competition for light between individuals, or stem structural stress due to snow load provide some plausible explanations.

Use of the Mechanistic Approach in Dendroecology
The application of our modelling approach to shrub-nitrogen data indicates that -as an alternative to traditional dendrochronological approaches -it may be useful to identify common mechanisms underlying the structure of 'multi-proxy' time-series data pertaining to a plant individual. Standard ring width methods involve the use of correlative approaches to identify a common signal across ring, isotope, or wood anatomy series, from many plants. Here we constructed an index of common mechanistic structure between individuals. As parameters in our models represent abstractions of plastic plant traits, a stand could be represented as a function of mechanistic limits to growth and envelopes of trait variability. Our approach allowed us to capture other controls to growth (e.g., temperature, precipitation) within estimated error terms, while focusing on N-limitation. We suggest that our approach may be expanded to: (a) include other Figure S2 Wood nitrogen concentration for ten Salix lanata L. shrub individuals.

Figure S3
Expected growth trajectories for all ten Salix lanata L. individuals given best-fitting parameters.
Methods S1 Ring width chronology development.

Methods S2 Creation of wood thin section 'cutting guides'.
Methods S3 Simulated annealing implementation.

Methods S4
Bristlecone F# scripts used for model-fitting and model-selection.
Notes S1 Sensitivity analysis of a global versus local calibration of the δ15N -N relationship.

Table S1
Parameter bounds from which starting values were drawn.

Figure 2
Schematic of time-series model structure that combines Salix lanata shrub ring-widths with stable nitrogen (N) isotope data to represent interactions between an individual shrub and soil N resources. This is the base model that is used to test all hypotheses. The filled boxes indicate measured time-series inputs, and hollow boxes the modelled stocks. Arrows and black labels indicate modelled N flows, alongside the processes that may give rise to them (blue text). Environmental losses (dotted line) may or may not be connected to soil N availability depending on the hypothesis under testing.

Figure 3
Dendroecological data for ten Salix lanata individuals, consisting of annual wood ring δ 15 N and annual ring increment relative to current radius. Standard error (red bound) was calculated as the standard deviation of replicate isotope samples (average of 0.303‰ over all δ 15 N series). Inset map shows spatial delineation between shrubs on the sloped valley side and lowland areas.

Accepted Article
This article is protected by copyright. All rights reserved 19.982 (both 3dp) fit to function in Section 3.3.4 using field measurements; B) stem volume and stem length; and C) stem basal radius and stem biomass. Relationships in B and C were generated using the allometric model from Niklas and Spatz (2004); the model was configured to create a shrubby form as specified in Supplementary Methods S4.

Data availability
The allometric, ring width, and stable isotope data that support the findings of this study are openly available at 10.5281/zenodo.4817326. In addition, the ring width and stable isotope time-series data have been submitted to the International Tree-Ring Data Bank and will be available at https://www.ncdc.noaa.gov/data-access/paleoclimatology-data/datasets/tree-ring.

Accepted Article
This article is protected by copyright. All rights reserved

Accepted Article
This article is protected by copyright. All rights reserved This article is protected by copyright. All rights reserved

Accepted Article
This article is protected by copyright. All rights reserved