Nitrogen restricts future treeline advance in the sub-arctic 1

Arctic environmental change has induced shifts in high latitude plant community composition and stature 9 with implications for Arctic carbon cycling and energy exchange. Two major components of high latitude ecosystems 10 undergoing change is the advancement of trees into treeless tundra and the increased abundance and size of shrubs. 11 How future changes in key climatic and environmental drivers will affect distributions of major ecosystem types is an 12 active area of research. Dynamic Vegetation Models (DVMs) offer a way to investigate multiple and interacting drivers 13 of vegetation distribution and ecosystem function. We employed the LPJ-GUESS DVM over a subarctic landscape in 14 northern Sweden, Torneträsk. Using a highly resolved climate dataset we downscaled CMIP5 climate data from three 15 Global Climate Models and two 21st century future scenarios (RCP2.6 and RCP8.5) to investigate future impacts of 16 climate change on these ecosystems. We also performed three model experiments where we factorially varied drivers 17 (climate, nitrogen deposition and [CO2]) to disentangle the effects of each on ecosystem properties and functions. We 18 found that treelines could advance by between 45 and 195 elevational meters in the landscape until the year 2100, de19 pending on the scenario. Temperature was a strong, but not the only, driver of vegetation change. Nitrogen availability 20 was identified as an important modulator of treeline advance. While increased CO2 fertilisation drove productivity in21 creases it did not result in any range shifts of trees. Treeline advance was realistically simulated without any tempera22 ture dependence on growth, but biomass was overestimated. As nitrogen was identified as an important modulator of 23 treeline advance, we support the idea that accurately representing plant-soil interactions in models will be key to future 24 predictions Arctic vegetation change. 25


Introduction 28
In recent decades, the Arctic has been observed to become greener (Epstein et  shows that 52% of treelines had advanced while the other half was stationary (47%), with only occasional instances of 48 retreat (1%) (Harsch et al., 2009). Similar patterns have been observed on the circumarctic scale, although latitudinal 49 treelines might be expected to shift more slowly than elevational treelines due to dispersal constraints (Rees et al., 50 2020). As trees close to the treeline often show ample storage of non-structural carbohydrates (Hoch and Körner, 2012) 51 it has been suggested that a minimum temperature requirement for wood formation, rather than productivity, constrains 52 treeline position (Körner, 2003(Körner, , 2015Körner et al., 2016). 53 Secondly, it has been hypothesised that indirect effects of warming might be equally or more important than direct 54 effects (Sullivan et al., 2015;Chapin, 1983). For example, rising temperatures and subsequently soil temperatures might 55 induce increased nitrogen mineralisation and plant nitrogen uptake (Chapin, 1983). Increased nitrogen uptake could in 56 turn enhance plant productivity and growth (Dusenge et al., 2019). Increased nitrogen uptake as a consequence of in-57 creased soil temperatures or nitrogen fertilisation have been shown to increase seedling winter survival among moun-58 tain birch (Betula pubescens ssp. tortuosa) seedlings (Weih and Karlsson, 1999; Karlsson and Weih, 1996). 59 Thirdly, experiments of elevated CO2 often show increased plant productivity and biomass increase, especially in trees 60 (Ainsworth and Long, 2005). Terrestrial biosphere models generally agree with this pattern (Hickler et al., 2008;Smith 61 et al., 2014;Piao et al., 2013). Although difficult to measure in field experiments, treeline position seems unresponsive 62 to increased [CO2] alone (Holtmeier and Broll, 2007). Whether treelines are responsive to increased productivity 63 through CO2 fertilisation might yield insights into whether treelines are limited by their productivity, i.e., photosynthe-64 sis, or ability to utilise assimilated carbon, i.e., wood formation. However, to what extent increased [CO2] drives long-65 term tree and shrub encroachment and growth remains poorly studied. 66 For treeline migration to occur, it is not only the growth and increased stature of established trees that is important, but 67 also the recruitment and survival of new individuals beyond the existing treeline (Holtmeier and Broll, 2007). Seedlings 68 of treeline species are sometimes observed above the treeline, especially in sheltered microhabitats (Hofgaard et al.,69 Gridcell vegetation and biogeophysical properties are calculated by averaging over a number of replicate patches, each 140 subject to the same climate forcing and intended to represent a random sample of the vegetation in the stand. Within 141 each patch, establishment, growth and mortality of individual tree or shrub cohorts are modelled annually (Smith et al. 142 2014). Establishment and mortality have both an abiotic (bioclimatic) and biotic (competition-mediated) component. 143 Vegetation dynamics, i.e. changes in the distribution and abundance of different PFTs in grid cells over time, are an 144 emergent outcome of the competition for resources between PFTs within an overall climate envelope determined by 145 bioclimatic limits for establishment and survival. Disturbance was accounted for by the occasional removal of all vege-146 tation within a patch with an annual probability of one per 300 years, representing random events such as storms, ava-147 lanches, insect outbreaks, and wind-throw. The study used three replicate patches within each 50x50m gridcell. 148 For summergreen PFTs we slightly modified the assumption of a fixed growing degree day (GDD) requirement for 149 establishment, using thawing degree days (TDD; degree days with a 0 °C basis; see Table S2.2) instead. 150

Historic period 152
A highly resolved (50x50m) temperature and radiation dataset using field measurements and a digital elevation model 153 (DEM) by Yang et al. (2011) provided climate input to the model simulations for the historic period . The 154 temperature in the lower parts of the Abisko valley in the dataset was influenced by the lake with milder winters and 155 less yearly variability as a result. At higher elevation, the temperature is more variable over the year and more depend-156 ent on the different solar angles between seasons (see Fig. S1.1; supplementary materials). 157 Monthly precipitation input was obtained from the Abisko Scientific Research Station weather records. Precipitation 158 was randomly distributed over each month using probabilities from the CRUNCEP v.7 dataset (Wei et al., 2014). We 159 assumed that local differences in precipitation can be neglected for our study domain and thus the raw station data was 160 used as input to LPJ-GUESS for the historic period. Nitrogen deposition data for the historic and future simulations 161 were extracted from gridcell including Abisko in the dataset produced by Lamarque et al. (2013). Nitrogen deposition 162 was assumed to be distributed equally over the study domain. 163 Data of soil texture was extracted from the WISE soil dataset (Batjes, 2005) for the Abisko area and assumed to be 164 uniform across the study domain. Callaghan et al (2013) reports that the soils around the Torneträsk areas are mainly 165 glaciofluvial till and sediments. Clay and silt fractions vary between 20-50% in the area (Josefsson, 1990) with higher 166 fractions of clay and silt in the birch forest and a larger sand content in the heaths. In the absence of spatial information 167 on particle size distributions, the soil was prescribed as a sandy loam soil with approximately 43% sand and approxi-168 mately equal fractions of silt and clay. 169

Future simulations 170
Future estimates of vegetation change were simulated for one low (RCP2.6) and one high (RCP8.5) emission scenario. 171 For each emission scenario, climate change projections from three global climate models (GCMs) that had contributed to the CMIP5 GCM ensemble (Taylor et al., 2012) were used to investigate climate effects on vegetation dynamics. The 173 GCMs (MIROC-ESM-CHEM, HadGEM2-AO, GFDL-ESM2M) were selected to represent the largest spread, i.e., 174 highest, lowest and near average, in modelled mean annual temperature for the reference period 2071-2100. Only mod-175 els that had contributed with simulations for both RCP2.6 and RCP8.5 were used in the selection. Monthly climate data 176 needed as input to LPJ-GUESS (temperature, total precipitation, and shortwave radiation) was extracted for the gridcell 177 including ANS for each GCM. 178 The historic climate dataset by Yang et al (2011) was extended into the projection period (2001-2100) using the delta 179 change approach, as follows. For each gridcell monthly differences were calculated between the projection climate and 180 the dataset by Yang et al. (2011) for the last 30-year reference period in our historic dataset . For tempera-181 ture, the difference was derived, while for precipitation and incoming shortwave radiation relative differences between 182 the two datasets were derived. The calculated monthly differences were then either added (temperature) to the GCM 183 outputs, or used to multiply (precipitation, radiation) the GCM outputs from 2001-2100, for each of the climate scenari-184 os used. Forcing data of atmospheric [CO2] for the two scenarios were collected from the CMIP5 project.  Table 1. 226

Determination of domains in the forest-tundra ecotone 227
In our analysis we distinguished between forest, treeline and shrub-tundra, defined as follows. Any gridcell containing 228 30% fractional projective cover or more of trees was classified as forest. This limit has been used by other studies in the 229 area (e.g., Van Bogaert et al., 2011) to determine the birch forest boundary. The treeline was then determined by first 230 selecting gridcells classified as forest. Any gridcell with 4 or more neighbours fulfilling the 30% cover condition crite-231 ria was classified as belonging to the forest. The perimeter of the forest was then determined through sorting out 232 gridcells with 4 or 5 neighbours classified as forest. Gridcells with fewer or more neighbors were regarded as tundra or 233 forest, respectively. Gridcells below the treeline were classified as forest in the analysis and gridcells above the treeline 234 were classified as tundra. 235 https://doi.org/10.5194/bg-2021-169 Preprint. Discussion started: 28 June 2021 c Author(s) 2021. CC BY 4.0 License.

Presentation of results 236
We present seasonal values for soil and air temperature. These are averages of DJF, MAM, JJA, and SON, referred to 237 as winter, spring, summer and autumn below. For the RCPs average values are presented with the ranges of the differ-238 ent scenarios within each RCP given in parenthesis. 239

Historic vegetation shifts 241
The dominating PFT in the forest and at the treeline was IBS which constituted 90% of the total LAI ( Fig. 2a-3a). The 242 only other tree PFT present in the forest was BINE, which comprised a minor fraction of total LAI. However, in the 243 lower (warmer) parts of the landscape BINE comprised up to 20% of total LAI in a few gridcells. Forest understory was 244 mixed but consisted mostly of tall and low evergreen shrubs and grasses. Shrub tundra vegetation above the treeline 245 was more mixed but LSE dominated with 51% of total LAI. Grasses comprised an additional 25% of total LAI and IBS 246 was present close to the treeline where it comprised up to 5% of LAI in some gridcells.  (Fig. 3b-c). The increase was more pronounced in RCP8.5, which also saw a large increase in low ever-269 green shrubs (LSE) at the end of the century (2090-2100). While the forest was still dominated by IBS, evergreen trees 270 (BNE and BINE) increased and together comprised approximately 15% of total LAI. The fraction of evergreen trees in 271 the forest correlated well with the degree of warming in each scenario. Forest productivity (GPP) was mainly driven by 272 tree PFTs and increased by 50% (12% -99%) for RCP2.6 and 177% (98% -270%) for RCP8.5. Above the treeline, low 273 shrubs (LSS and LSE) contributed most to annual GPP increases, which increased by 33% (-12% -67%) and 239% 274 (105% -372%) in RCP2.6 and RCP8.5, respectively. The productivity increase translated into a biomass C increase of 275 the same magnitude in both the forest and above the treeline. 276 The average summer air temperature at the treeline between the last decade in the historic and projection periods in-277 creased by 0.3°C and 6.7°C for the coldest (GFDL-ESM2M-RCP2.6) and warmest (MIROC-ESM2M-RCP8.5) GCM 278 scenario, respectively. The advance of the 6°C JJA soil temperature isoline was more rapid than the treeline advance 279

Model experiments 284
A slight treeline advance at the end of the projection period (2090-2100) of approximately 11 elevational meters was 285 seen in the control simulation where all drivers were held constant. This revealed a lag from the historical period, likely 286 resulting from smaller trees that had established in the historic period that matured during the projection period. 287

Climate change 288
Treeline advance occurred in all climate change scenarios although the rate was not uniform throughout the projection 289 period (Fig. 5). When driven by climate change alone, migration rates were faster than the projection simulations where 290 nitrogen deposition and CO2 were also changed (section 3.2). Treeline advance in climate change only scenarios ranged 291 between 60 elevational meters (HadGEM2-AO-RCP2.6) to 245 elevational meters (MIROC-ESM-CHEM-RCP8.5) 292 over the 100 year projection period. 293 Tree productivity was strongly enhanced by air temperature over the whole study domain (Fig. 6a). Other climate fac-294 tors such as precipitation and net shortwave radiation also correlated with productivity, however, these correlations 295 were weaker (Fig. S1.5; S1.6; supplementary materials). Annual precipitation increased in all climate change scenarios 296 (Table 2). In the lower parts of the valley, the increased precipitation did not result in any increased soil moisture during 297 summer as losses through evapotranspiration driven by temperature exceeded the additional input. Spring and autumn 298 soil moisture increased in the forest, mainly because of earlier snowmelt and thawing ground in spring and a relatively 299 weaker evapotranspiration in autumn. Above the treeline, soil moisture increased as the lower temperatures and LAI did 300 not drive evapotranspiration as strongly as in the lower parts of the valley and the increased rain thus outweighed the 301 slightly increased evapotranspiration. Increased tree productivity in the forest resulted in an increased LAI of 18-90%, which in absolute terms was equivalent 303 to an increased LAI of 0.3-1.5 m 2 m -2 . BNE appeared in the forest and dominated in a few gridcells. In most places 304 BNE constituted approximately 5% of total LAI. Tall shrub (HSE and HSS) productivity and LAI increased in the for-305 est, however, the increase was negatively correlated with temperature, i.e., increase was highest in the coldest climate 306 change scenarios. Above the treeline, tall shrubs showed the opposite pattern, increasing by 8-50% to finally constitute 307 10-36% of total LAI. 308 Higher soil moisture content in spring and autumn favoured trees in the whole ecotone, while forest understory suffered 309 from earlier onset of the growing season with subsequent flushing of the leaf and light shading from taller competitors. Radiation correlated positively with the growth of tree PFTs, with spring and autumn radiation found to be especially 321 important for height and biomass increase (Fig. S1.7; supplementary materials). Increased radiation provided a competi-322 tive advantage for taller trees and shrubs to shade out lower shrubs and grasses in the forest. Shrubs above the treeline 323 were also favoured by increased light. 324 Net nitrogen mineralisation, i.e., the difference between microbial nitrogen mineralisation and immobilisation, at the 325 treeline showed great variation between different climate change scenarios, ranging from a 4% decrease in one scenario 326 (GFDL-ESM2M-RCP8.5) to a 79% increase in the strongest warming scenario (MIROC-ESM-CHEM-RCP8.5). In 327 absolute terms, the latter increase corresponds to the nitrogen load in the 7.5x increased nitrogen deposition scenario. 328 Interestingly, despite very different plant available nitrogen and warming, the two scenarios displayed a similar result-329 ing (2090-2100) treeline elevation (Fig. 5a). 330

CO2 331
Productivity increase in most PFTs correlated well with a [CO2] increase (Fig. 6b). Total GPP averaged over the forest 332 increased between 2-10% depending on the [CO2] scenario, with the largest increase in RCP8.5 and smallest in RCP2.6. 333 The CO2 fertilisation effect was not uniform within the landscape, but stronger towards the forest edge with increases 334 from 2% to 18% from the weakest to the strongest [CO2] scenario. IBS NPP increased uniformly over the forest with 335 2.5-8.4% but decreased above the treeline. Thus, the productivity of the two dominant PFTs (IBS in the forest and LSE 336 above the treeline) was reinforced in their respective domains. The increased productivity translated into a 1-5% in-337 https://doi.org/10.5194/bg-2021-169 Preprint. Discussion started: 28 June 2021 c Author(s) 2021. CC BY 4.0 License. crease in tree LAI in the forest while low shrub LAI increased with 24-77%. Likewise, increase in leaf area of low 338 shrubs was largest on the tundra under elevated [CO2], which saw a 15-40% LAI increase in the low and high [CO2] 339 scenario respectively. Above the treeline, the productivity of grasses and low shrubs responded strongly to the CO2 340 fertilisation with a 350% increase in GPP for grasses and 150% increase for low shrubs. The additional litter fall pro-341 duced by the increased leaf mass did not lead to any increase in N mineralisation. However, immobilisation of nitrogen 342 through increased uptake by microbes increased with 2-6% between the lowest and highest [CO2] scenarios, yielding a 343 net reduction of plant available nitrogen. Furthermore, the productivity increase did not drive any range shift of the 344 forest, i.e., the treeline remained stationary in all [CO2] scenarios (Fig. 5b). 345

Nitrogen deposition 346
Productivity of woody PFTs was in general positively correlated with nitrogen concentration in the different nitrogen 347 deposition scenarios. In contrast, productivity of grasses was negatively correlated (Fig. 6c) as they suffered in the light 348 competition with the trees. Annual GPP of trees (especially IBS) was positively correlated throughout the whole eco-349 tone, however, the increase in GPP was larger towards the forest boundaries than in the lower parts the forest when 350 nitrogen was added. Nitrogen stressed plants in the model allocate more carbon to their roots at the expense of foliar 351 cover when they suffer a productivity reduction. In the two scenarios with decreasing nitrogen deposition (RCP2.6; 352 RCP8.5) there was an overall reduction in LAI in both the tundra and the forest of a magnitude 6-10%. The largest 353 reduction was seen in tree PFTs, which have the largest biomass and consequently could be assumed to have the highest 354 nitrogen demand, followed by tall shrubs. Low shrubs and grasses did however increase their LAI in the forest when 355 nitrogen input decreased resulting from a decreased light competition from trees. Above the treeline, LAI of low shrubs 356 and grasses PFTs also decreased with less nitrogen input. 357 In all scenarios with increasing nitrogen deposition there was an advancement of the treeline in the order of 10-85 ele-358 vational meters with smallest (2x nitrogen deposition) having the smallest change in treeline elevation and vice versa 359 for largest input (10x nitrogen deposition) scenario (Fig. 5c). In the scenarios where nitrogen input was constant or 360 decreasing, the treeline remained stationary. 361

Discussion 362
In our simulations, rates of treeline advance were faster under climate change-only scenarios than when all drivers were 363 changing. This revealed nitrogen as a modulating environmental variable, as nitrogen deposition was prescribed to 364 decrease in both the RCP2.6 and RCP8.5 scenarios. In contrast to previous modelling studies of treeline advance (e.g., 365 Paulsen and Körner, 2014), we include not only temperature dependence on vegetation change, but also the full nitro-366 gen cycle and CO2 fertilisation effects . Increased nitrogen deposition induced treeline advance, 367 further illustrating the importance of nitrogen dynamics in our results. In the elevated [CO2] scenarios, higher produc-368 tivity in all plants was induced, but productivity enhancements alone did not lead to significant treeline advance. How-369 ever, enhancement of productivity in combination with an allocation shift to the plant canopies, enabled by a greater 370 nitrogen uptake, favoured taller plants over their shorter neighbours in the competition for light within the model. Field 371 experiments with nitrogen fertilisation have shown that mountain birch at the treeline displays additional growth after nitrogen additions (Sveinbjörnsson et al., 1992). Furthermore, fertilisation with nitrogen improved birch seedling sur-373 vival above the treeline (Grau et al., 2012), and is thus likely important for establishment and growth of new individuals 374 to form a new treeline. This pattern is in line with the hypothesis that productivity increase alone does not drive range 375 shifts of trees (Körner, 2015). As has also been pointed out by others (Hofgaard et  shrub productivity and growth consequently are responsive to CO2 fertilisation. 450

Conclusions 451
In this study we identified nitrogen cycling and availability as an important modulator of treeline advance. Internal 452 cycling of nutrients in soils is the main source of nutrients for Arctic plants (Chapin, 1983). The model performed well 453 regarding rates of shrub increase and treeline advance but overestimated biomass carbon in the forest. Treeline migra-454 tion rates were realistically simulated although the model did not represent temperature limitations on tree growth. 455 While a decoupling between productivity and growth in the model could potentially have improved estimates of bio-456 mass carbon, it was not needed to correctly predict treeline elevation. Instead, our results point to the importance of 457 indirect effects of rising temperatures on tree range shifts, especially with regards to nutrient cycling. Furthermore, soil 458 moisture strongly influenced vegetation composition within the model with implications for treeline advance. How 459 models represent nutrient uptake and cycling, as well as a better empirical understanding of processes that determine 460 tree and shrub growth will be key to better predictions of Arctic vegetation change and carbon and nitrogen cycling. 461 Models are a valuable aid in judging the relevance of these processes on the pan-Arctic scale.