Increasing tree density accelerates stand‐level nitrogen cycling at the taiga–tundra ecotone in northeastern Siberia

Abstract

As climate warms, tree density at the taiga–tundra ecotone (TTE) is expected to increase, which may intensify competition for belowground resources in this nitrogen (N)‐limited environment. To determine the impacts of increased tree density on N cycling and productivity, we examined edaphic properties indicative of soil N availability along with aboveground and belowground tree‐level traits and stand characteristics related to carbon (C) and N cycling across a tree density gradient of monodominant larch (Larix cajanderi) at the TTE in far northeastern Siberia. We found no consistent evidence from soil, tree, or stand‐level N cycling characteristics of lower N availability or greater intraspecific competition for N with increased density. Active layer thickness declined, but resin‐sorbed N and soil organic layer thickness did not covary with increased tree density. There was, however, greater allocation belowground to stand‐level coarse and fine roots with increased tree density, an allocation pattern suggestive of limited soil resources. Foliar traits related to C (%C, δ¹³C, and resorption) were responsive to density indicating the importance of non‐nutrient resources, like light, to foliar stoichiometry. As tree density increased and individual trees had lower productivity, tree‐level N and biomass pools aboveground and belowground declined tracking decreases in N uptake, N resorption, N use efficiency, and allocation to slow cycling tissues like wood. At the stand level, our findings show high N turnover with increased N acquisition, allocation to short‐lived tissues with relatively high N content and reduced N residence time, and greater stand productivity as tree density increased. Yet, these positive relationships were curtailed at the highest tree densities. Our observations of shifts in biomass, C and N allocation, and loss aboveground, along with greater root density with increased tree density, could have strong impacts on C and N cycling and should be represented in models of TTE dynamics and feedbacks to climate.

INTRODUCTION

As high-latitude ecosystems warm, arctic forests at the taiga–tundra ecotone (TTE; Montesano et al., 2020) are expected to change in structure and expand into the treeless tundra (Grace et al., 2002; Shevtsova et al., 2020), thus impacting the climate system through increased ecosystem carbon (C) storage and decreased albedo (Chapin III et al., 2005; Pearson et al., 2013). The TTE is a transition zone from forested to tundra landscapes at the cold edge of the boreal forest and encompasses gradients in canopy cover, density, height, and aboveground biomass (Montesano et al., 2020). Increased tree density at the TTE can occur due to the direct effects of warming (Holtmeier & Broll, 2007; MacDonald et al., 2008) or indirect effects of warming on disturbance regimes like wildfire that impact tree establishment and growth (Brown, 2010; Landhausser & Wein, 1993). Historically, changes in TTE dynamics have been attributed to the warming of growing season (Korner & Paulsen, 2004) or winter air temperature (Harsch et al., 2009). Contrary to expectation, only half of the global TTE sites show forest advance with warming (Harsch et al., 2009; Rees et al., 2020), indicating that other mechanisms like nutrient availability (Ellison et al., 2019; Gruber et al., 2018; McNown & Sullivan, 2013; Sullivan et al., 2015) could be important to tree growth and potential biome shifts. The formation of high-density TTE forests, a precursor to forest advance, is likely to increase competition for limited resources (Liang et al., 2016; Wang et al., 2016) like soil nutrients that impact boreal tree growth (Högberg et al., 2017). Therefore, to accurately predict TTE dynamics as the climate warms, we need to understand the effect of increased density on tree- and stand-level acquisition and use of limited belowground resources and its coupling to productivity.

Soil nitrogen (N) supply constrains plant growth in the boreal forest (Högberg et al., 2017; Tamm, 1991) and tundra (Mack et al., 2004). The cold soils of the TTE are underlain by permafrost (Brown et al., 1997). Trees, therefore, have limited soil volume from which to acquire nutrients, and soil nutrient availability is low due to temperature constraints on decomposition (Schuur & Mack, 2018). Warming across the Arctic may ameliorate N limitation due to increased N mineralization in active layer soils, the seasonally thawed portion of the soil profile, and N release from deeper, thawed, permafrost soils (Salmon et al., 2018). However, N supply may vary with tree density because of biophysical feedbacks between vegetation cover, soil thermal regimes, and nutrient cycling (Bonan, 1990; Kropp et al., 2019). As tree density increases, soils become increasingly shaded, reducing surface soil temperature, active layer thickness, decomposition rates, and thus N supply.

Increased tree density could represent a scenario of low soil N supply, high N demand, and thus high intraspecific competition (Wang et al., 2016; Wieczorek et al., 2017) resulting in nutrient limitation of growth. In low-N environments, plants may have lower foliar %N, more depleted foliar δ¹³C and δ¹⁵N signatures (Craine et al., 2009; Li et al., 2007), and low wood and leaf production (Dybzinski et al., 2011; Litton et al., 2007). They may also have high N use efficiency (NUE) due to high N productivity (g biomass g⁻¹ N) or longer N mean residence time (NMRT; Chapin III, 1980, Lambers & Poorter, 1992), and increased allocation to N uptake, including allocation to root production or mycorrhizal fungi as signaled by greater root %N and δ¹⁵N (Hobbie & Colpaert, 2003; Hobbie & Hobbie, 2008; Litton et al., 2007). Examination of soil conditions, individual tree traits, and stand-level characteristics together can therefore inform our understanding of the relationship between increased tree density, N cycling, productivity, and the potential for nutrient constraints on TTE dynamics.

To determine how tree density impacts forest C and N cycling at the TTE, we measured aboveground and belowground metrics related to N cycling and productivity along a gradient of increasing larch (Larix cajanderi) tree density in northeastern Siberia. The gradient analysis approach allows not only for both the detection of change in ecological structure and function in association with the environment, but also for the extrapolation of observations across multiple spatial scales (Gosz, 1992), that is, from trees to stands to biomes. Our findings, therefore, can inform conceptual models of TTE dynamics with climate warming. We asked how tree density impacts:

1. edaphic factors indicative of soil N availability, including resin-sorbed soil N in organic and mineral horizons, and active layer and soil organic layer (SOL) thickness,
2. stand-level larch coarse-root density, proportion of fine roots, root tissue C and N content and isotopes, and tree-level root biomass, which provide metrics of allocation to the acquisition of soil resources and N availability, and
3. aboveground tree-level traits and stand-level characteristics related to N uptake, allocation, and loss.

METHODS Site description

We conducted our research in the larch forests around the Northeast Science Station (NESS; 68.74° N, 161.40° E) in northeastern Siberia close to Cherskiy, Sakha Republic, Russian Federation. The NESS is located on the Kolyma River, ~250 km north of the Arctic Circle and ~130 km south of the Arctic Ocean. Climate is continental, with warm summers (June average = 12°C), cold winters (January average = −33°C), and average annual temperature of −11.6°С (Cherskiy Meteorological Station; https://rp5.ru/Weather_archive_in_Cherskiy). Annual precipitation is low (210 mm year⁻¹; Zimov, unpublished data), with approximately half falling during summer (Berner et al., 2013). Forests in this region of the Russian Far East are typically open-canopied, sparse stands dominated by L. cajanderi (Alexander et al., 2012, 2018; Berner et al., 2012), a deciduous needleleaf conifer, which grows in habitats with continuous permafrost and a short growing season (Abaimov, 2010). Trees in this region are short in stature, generally <10-m tall, and forests have relatively low aboveground biomass (<1200 g m⁻²) compared to those further south (Berner et al., 2012). Stands of larch in our study system are underlain by carbon-rich, Pleistocene-age Yedoma permafrost with high ice content (Strauss et al., 2017). The active layer thickness across the stands is impacted not only by climate but is also sensitive to tree density and canopy cover more so than variation in SOL thickness unlike the ecosystem-protected permafrost (Shur & Jorgenson, 2007) further south. While soil temperature declines, soil moisture remains similar as thaw depth declines and tree density increases (Loranty & Alexander, 2021; Appendix S1: Figures S1, S2). Understory vegetation is composed of tall (>1 m) deciduous and short (<1 m) mostly ericaceous shrubs, forbs, grasses, lichens, and mosses. Across the stands we studied, tall shrubs declined in abundance, while mosses increased with greater tree density (Paulson et al., 2021).

Field sampling

To determine whether C and N cycling metrics at the individual tree or stand level vary with tree density, we sampled 26 larch stands across a tree density gradient from 2010 to 2017. The tree productivity, understory composition, ecohydrology, SOL depth, and active layer depth for many of these sites have been described previously (Kropp et al., 2019; Paulson et al., 2021; Walker et al., 2021). Tree densities in these sites ranged from 0.03 to 3.70 trees m⁻². Percent canopy cover (Paulson et al., 2021) and tree basal area (3–19 m² ha⁻¹) increased with tree density, while average tree diameters were greater in low-density (5 cm) compared to high-density (3 cm) stands (Walker et al., 2021). All stands were within the perimeter of a 1940 fire scar located ~2 km from the NESS, and thus, trees were of a similar age (Walker et al., 2021). The gradient of tree densities has been attributed to the impact of soil burn severity on larch recruitment, where high fire severity coupled with adequate seed source resulted in high recruitment density and subsequently high tree density at maturity (Alexander et al., 2018; Paulson et al., 2021). Variation in canopy cover in the gradient (Paulson et al., 2021) mirrors that observed across the TTE (Montesano et al., 2020). We sampled soil and vegetation parameters (Appendix S1: Figure S3) in three plots located at least 30 m apart within each of the 26 stands (each stand ~0.5 ha). Plots consisted of a variable-width, 30-m-length belt transect. The width of the belt transect ranged from one m wide in the stands with the highest tree density to eight m wide in the stands with the lowest tree density (Paulson et al., 2021; Walker et al., 2021). Across the gradient, climate, regional species pool, topography, parent material, and stand age are relatively constant, allowing us to test the effects of tree density on soil metrics, traits of individual larch trees, and characteristics at the stand level related to C and N cycling.

Active layer and soil organic layer

In each stand, we measured the September thaw depth in 2017, a close proxy of active layer thickness, at five points running down the center of each plot (Appendix S1: Figure S3). We measured the SOL thickness in July–August 2015 in each stand at two locations (~30 m apart) adjacent to each of the three sampling plots (n = 6 soil samples/stand; Appendix S1: Figure S3).

Available soil nitrogen

We estimated NO₃⁻ and NH₄⁺ availability in the field using mixed-bed ion exchange resins (J.T. Baker Scientific). We installed one bag with 10-g moist weight resin in the organic layer and another in the upper mineral soil (~5 cm below the surface of each soil layer) at two locations ~30 m apart adjacent to two sampling plots (n = 2 resin bags horizon⁻¹ plot⁻¹ × 2 plots stand⁻¹; Appendix S1: Figure S3). Resin bags were incubated for 1 year (July 2015–July 2016). Resins were extracted in 50 ml 2 M KCl by shaking for 4 h, filtered through a glass fiber filter (grade GF/A), measured on a SmartChem 200 (Unity Scientific, Brookfield, CT, USA), and expressed in milligrams of N per year. Only a subset of extracts were successfully exported due to travel restrictions, and thus, we have resin-sorbed N measurements from 9 of 26 stands.

Stand root characteristics

We quantified coarse (≥2 mm in diameter) and fine (<2 mm in diameter) root density (in grams per cubic centimeter) in the SOL within each stand from excavated soils collected in July and August 2015–2017 at two locations ~30 m apart adjacent to two of the sampling plots in each stand (n = 4 root samples/stand, n = 25 stands; Appendix S1: Figure S3). All sampling points were at least 30 cm away from a mature tree, but distance to the nearest tree was highly variable across stands because of the variation in tree density, with many more nearby trees in high- and moderate-density stands compared with low-density stands. Coarse roots were excavated from the SOL using a soil saw within a 0.04-m² area to ~20 cm depth. Fine roots were obtained from a subsample (~10 cm length × 10 cm width) of the coarse root soil sample. Organic soil was removed from coarse and fine roots by gently washing under running water. Live roots were determined by tensile strength and color. Coarse live roots were sorted by type (larch or other), while root species was not distinguished for fine roots. Roots were dried at 60°C for 48 h, and stand-level root densities were calculated by dividing the dry mass of the roots by the soil volume sampled (Alexander et al., 2021b). As a proxy for L. cajanderi fine-root density, we estimated fine-root density as half the coarse larch root density in the organic horizon based on estimates of a congeneric L. gmelnii growing on permafrost soil in Siberia (Kajimoto et al., 1999). We calculated the larch proportion of total fine-root density as the estimate of larch fine-root density divided by the total fine-root density. To assess tree-level allocation to fine roots, the estimate of stand-level fine-root density was converted to biomass in grams per square meter and divided by tree density (in trees per square meter).

To determine the C and N content and isotope signatures of larch fine roots, we sampled larch fine roots from 10 of the 26 stands that were representative of the range of densities (10 stands 0.10–3.70 trees m⁻² and 26 stands 0.03–3.70 trees m⁻²) across the gradient in September 2017. Root δ¹³C and δ¹⁵N can reflect C assimilates delivered to the roots, fungal colonization or biomass, or N source (Brüggemann et al., 2011; Hobbie & Colpaert, 2003; Templer et al., 2007) and thus are suggestive of the resource economics of the stand. We collected roots in each stand from organic and mineral soils at six sampling locations located ~10 m apart from one another and different from those described above to obtain root densities. At each sampling location, we removed a 10 × 10 cm organic soil monolith with a depth of the full organic horizon. Below each monolith, we sampled the top 20 cm of mineral active layer soils at 10-cm increments. From each stand, we pooled the six soil samples by horizon and depth increment and processed as one sample. We sampled larch roots from each pooled sample. Larch fine roots were not found in appreciable amounts (<0.003 g cm⁻³) in mineral soils; however, we analyzed the root chemistry of the root fragments recovered. Roots were dried at 60°C for 48 h, ground, and run for C and N concentration and isotope analysis on an elemental analysis-isotope ratio mass spectrometer (Delta Advantage, Thermo Fisher Scientific, Waltham, MA, USA) coupled to an elemental combustion analyzer (Costech ECS4010, Valencia, CA, USA).

Aboveground tree and stand characteristics

We assessed 16 foliar traits of individual trees based on C and N concentrations and δ¹³C and δ¹⁵N isotopes of foliage sampled in July and September in relation to tree density (Appendix S2: Table S1). In July 2017 during peak biomass, we sampled foliage from five overstory larch trees within each of the 26 stands in the density gradient (Appendix S1: Figure S3). We selected a tree near the 0-m mark of each plot and a tree near the 30 m for two of the plots, ensuring that all trees were at least 30 m apart. The trees represent the size distribution of trees in the stand as a whole. In September 2017, we resampled foliage at senescence that had yet to fall as litter from the same five trees. Leaves were collected from lower branches within reach of the ground, dried at 60°C for 48 h, ground, and run for C and N concentration and isotope analysis as described above for roots. The foliar C and N concentrations were also utilized in calculations of other tree traits and stand metrics (Appendix S2: Table S1).

In addition to foliar traits of individual trees, we characterized other aboveground C and N cycling metrics at the individual tree level and the stand level. In total, we calculated 49 tree-level and 24 stand-level characteristics (i.e., C, N, and biomass pools, resorption, N uptake, N production, NMRT, and NUE, calculations shown in Appendix S2: Table S1). We used inventory data collected from 2010 to 2017 at the three plots located within each of the 26 stands (Alexander et al., 2021a). In brief, we measured diameter at breast height (≥1.4-m tall) or basal diameter (<1.4-m tall) for each live larch tree within each plot (i.e., belt transect with larger area for lower density). Estimates of larch aboveground biomass were based on allometric equations (Alexander et al., 2012), and production was based on the 10-year average ring width measurements obtained from basal cores or disks ~30 cm above the forest floor from 5 to 10 trees per stand (Walker et al., 2021). Our calculations of C and N cycling metrics presented in this study are based on these biomass and productivity values for individual trees and measurements of C and N content of tree tissues. We calculated stand-level metrics by summing the tree-level pool estimates for each stand (Appendix S2: Table S1).

Statistical analysis

All statistical analyses were conducted in R 3.3.2 (R Core Team, 2019).

Relationships between tree density, soil parameters, and N availability

We tested whether active layer thickness was related to the fixed factor tree density using a linear mixed effects model with stand and plot nested in stand as random factors with the package “lme4” (Bates et al., 2015). We tested whether the thickness of the soil organic horizon varied with the fixed factor tree density with stand and plot nested in stand as random factors. In this and all subsequent linear mixed effects models, the random factor was used to induce a correlation structure between observations from the same stand or plot within stand (Zuur et al., 2009) and allows different intercepts for each stand or plot. Response variables were transformed (log or square root) when necessary to meet model assumptions of normality and homogeneity of variance, which were examined with diagnostic and residual plots. The significance of fixed effects was determined using maximum-likelihood ratio tests comparing the full model to the null model from which we report p values and associated test statistics and confirmed using Akaike information criterion (ΔAIC > 2.0; Burnham & Anderson, 2002; Zuur et al., 2009). For all models with more than one fixed factor, we tested the significance of interactions and then reduced the model to an additive structure if possible.

To determine whether nitrogen availability varied with tree density, we tested whether NO₃⁻-N, NH₄⁺-N, NH₄⁺-N:NO₃⁻-N, or dissolved inorganic nitrogen (DIN; NO₃⁻-N + NH₄⁺-N), all indicators of N availability and cycling rates, related to the fixed factor tree density and soil horizon (organic or mineral) using separate linear mixed effects models for each N form with stand and plot nested in stand as random factors. For these and all subsequent models where we tested multiple response variables with the same experimental units in relation to the fixed factors, we corrected p values by controlling for false discovery rates with the Benjamini–Hochberg procedure (Benjamini & Hochberg, 1995). We tested for significant differences between organic and mineral horizons with the post hoc test with the Tukey adjustment in the “emmeans” package (Lenth et al., 2019).

Relationships between tree density and stand-level root characteristics

We constrained our analysis of larch root density to the organic horizon where coarse and fine roots were most abundant. Coarse-root density was modeled in relationship to the fixed factor tree density with stand and plot nested in stand as random factors. Two high-density coarse root samples were removed from statistical analysis due to a disproportionate influence on the data distributions and model outcomes (Appendix S3: Figure S1). The proportion of larch fine-root density was modeled in relation to the fixed factor tree density with stand and plot nested in stand as random factors. In nine out of 100 samples, our literature-based estimate of fine-root density as half of coarse-root density yielded a proportional value over one, and these samples were removed from analysis. The standardized fine root biomass (in grams of fine root biomass per tree) was tested in relation to the fixed factor tree density with stand and plot nested in stand as random factors. We tested whether fine root δ¹³C, δ¹⁵N, %C, %N, and C:N varied with tree density and soil horizon using a linear model. We compared differences between soil horizons with a post hoc test using the Tukey adjustment described above.

Relationships between tree density and aboveground tree and stand characteristics

We quantified how larch foliar traits varied with the fixed factor tree density with stand as the random factor. Due to the observed late-season relationship between foliar δ¹⁵N and tree density, we further explored whether foliar δ¹⁵N from needles collected at senescence in September varied with the fixed factors tree density and percent N resorbed with stand as a random factor. We assessed whether the relationship between foliar δ¹⁵N and δ ¹³C was similar at peak biomass in July and at senescence in September using two models. First, we tested whether the response variable foliar δ¹⁵N covaried with the fixed factors δ¹³C in July and tree density and again tested this with the foliar signatures from foliage collected in September. In both models, we accounted for stand as a random factor. We also calculated Pearson correlation coefficients between foliar δ¹⁵N and δ ¹³C in July and again for September to determine whether the relationship between C and N isotope signatures was consistent across seasons.

We assessed whether 49 tree-level C and N cycling metrics (Appendix S2: Table S1), beyond the foliar traits described above, varied in response to the fixed effect tree density with stand and plot nested in stand as random factors. We conducted the same analyses at the stand level for 24 characteristics (Appendix S2: Table S1) using linear regression with each stand-level C and N cycling metric as a response variable and tree density as the predictor variable. Graphical analysis of tree and stand aboveground metrics in relation to density suggested a curvilinear response; we, therefore, tested whether a quadratic function improved linear model fit. We evaluated the AIC between the linear and the quadratic models. When the linear model AIC was greater than the quadratic model AIC by two and the orthogonal quadratic parameter estimates indicated additional significance to the linear term, we report the quadratic model. At the tree and stand levels, we graphically present seven metrics commonly used to assess N cycling (uptake, allocation, and loss): N pool in biomass (stems, branches, and foliage) and aboveground net primary production (ANPP) (woody ANPP and foliage), N uptake (ANPP N–N in September foliage), N resorbed (peak biomass July foliage N pool − senescent September foliar N pool), NMRT (N pool in biomass/N uptake), N productivity (ANPP biomass pool/ANPP N pool), and NUE (NMRT × N productivity).

RESULTS Relationships between tree density, soil parameters, and N availability

Across the 26 stands, active layer thickness decreased with increasing tree density (F = 5.62, p = 0.02, R²_m = 0.11, and R²_c = 0.58; Figure 1). The minimum measured active layer thickness was 20 cm, and the maximum was 122 cm. In contrast, SOL thickness was highly variable (3–21 cm; mean: 11 cm) but did not significantly vary with tree density (F = 0.61, p = 0.44).

Generate image description

FIGURE 1. Active layer thickness across 26 stands with increasing tree density near the Northeast Science Station, Cherskiy, Russia. Points are raw data, and the line is the linear mixed effects model fit with the 95% confidence interval shaded in gray. Darkness of green scales from low to high density.

Overall, we found that resin-sorbed N did not vary significantly with tree density. We observed that forms of resin-sorbed N varied from 0 to 15.04 mg N year⁻¹ in organic and mineral soils (Appendix S2: Table S2). Nitrate was related to soil horizon (F = 7.84, p = 0.02) but not tree density (F = 0.00, p = 1.00) and was higher in the mineral horizon than in the organic horizon (t = 2.80, p < 0.01). Neither NH₄⁺, NH4⁺-N:NO3_--N, nor DIN were significantly related to soil horizon or tree density (p > 0.05; Appendix S2: Table S2).

Relationships between tree density and stand-level root characteristics

Larch coarse-root density ranged from 0 to 0.017 g cm⁻³ and increased with tree density (F = 5.87, p = 0.02, R²_m = 0.09, and R²_c = 0.33; Figure 2a). The proportion of fine-root density attributed to larch roots ranged from 0.00 to 0.93 of the total root density, and fine-root density ranged from 0 to 0.007 g cm⁻³. The proportion of fine-root density attributed to larch roots increased with greater tree density (F = 25.94, p < 0.01, R²_m = 0.28, and R²_c = 0.37; Figure 2b). The tree-level estimate of fine root biomass declined with increased tree density (F = 13.29, p < 0.01, R²_m = 0.18, and R²_c = 0.33; Figure 2c). Of note, exclusion of the two highest biomass samples did not influence model interpretation (F = 19.76, p < 0.05, R²_m = 0.18, and R²_c = 0.18).

Generate image description

FIGURE 2. Larch coarse root density (a), proportion of total fine-root density (b), and fine root biomass (c) in relation to increasing tree density. Panels (a) and (b) depict stand-level estimates, while panel (c) depicts tree-level estimates. Points are raw data, and the line is the linear mixed effects model fit with the 95% confidence interval shaded in gray. Inclusion of the two highest biomass samples in panel (c) did not affect the interpretation of model significance. Darkness of green scales from low to high density.

Larch fine root C:N did not vary with tree density (F = 0.13, p = 0.72) but did vary by soil horizon (F = 6.00, p = 0.02). Roots had higher %C (est = 6.76, t ratio = 5.02, p < 0.001) and %N (est = 0.23, t ratio = 5.75, p < 0.001) in organic compared with mineral horizons, but C:N was lower in the organic than in the mineral horizon (t ratio = −2.41, p = 0.03). The δ¹³C isotope signature of fine roots did not vary with tree density or by soil horizon (F = 0.63, p = 0.54), while the δ¹⁵N signature varied by horizon (est = −1.51, t ratio = −2.54, p = 0.02), with roots from the organic horizon showing less δ¹⁵N enrichment (organic 0.47 ± 0.44‰ and mineral 1.98 ± 0.40‰).

Relationships between tree density and aboveground tree and stand characteristics

Foliar C traits of individual trees were related to variation in tree density, while most foliar N traits were not (Figure 3a–f; Appendix S2: Table S3). Foliage had greater %C in high-density stands in July during peak biomass, but lower %C in high-density stands by September at the time of senescence (Figure 3a,b; Appendix S2: Table S3). Foliar C resorption increased with tree density (Figure 3c), while, on average, larch resorbed 73.43 ± 0.84 SE %N, but this did not vary with tree density. The foliar δ¹³C in July (Figure 3d) and the seasonal change in δ¹³C (July–September; Figure 3e) declined with increasing tree density. However, the foliar δ¹³C signatures in low-density stands were more depleted than those in high-density stands by the end of the growing season in September compared to July (Appendix S3: Figure S2).

Generate image description

FIGURE 3. Foliar traits in relation to increasing tree density. Points are raw data, and the line is the linear mixed effects model fit with the 95% confidence interval shaded in gray. Inclusion of the highly depleted δ15N sample in panel (f) did not affect the interpretation of model significance. Darkness of green scales from low to high density.

During the peak of the growing season, foliar N isotope signatures of individual trees were not related to tree density (Appendix S2: Table S3). However, by September, the foliar δ¹⁵N signature decreased as tree density increased (Figure 3f; Appendix S2: Table S3). As a greater percentage of N was resorbed, we observed greater isotopic depletion in the September leaves (F = 4.44, p = 0.04, R²_m = 0.17, and R²_c = 0.52; Figure 4a). Trees with greater foliar δ¹⁵N enrichment also showed enrichment in δ¹³C in July (F = 5.66, p = 0.02, R²_m = 0.06, and R²_c = 0.49; Figure 4b) and September (F = 3.61, p = 0.06, R²_m = 0.17, and R²_c = 0.52; Figure 4c) while accounting for tree density. September δ¹³C and δ¹⁵N signatures showed high fidelity with those observed in July (δ¹³C Pearson's r = 0.52; δ¹⁵N Pearson's r = 0.57).

Generate image description

FIGURE 4. Relationships between larch foliar traits across the tree density gradient. Points are raw data, and the line is the linear mixed effects model fit with the 95% confidence interval shaded in gray. The foliar δ15N signatures in September declined with greater N resorbed (a) and increased with greater foliar δ13C enrichment at peak biomass in July (b) and at the time of senescence in September (c). Inclusion of the highly depleted September foliage sample did not affect the interpretation of model significance. Darkness of green scales from low to high density.

In general, metrics of C and N cycling measured on the per tree basis declined with tree density (32 of 49 traits; Appendix S2: Table S3; Figure 5a–h), including N pools, N uptake, N resorption, and NUE (Figure 5a–d,g). There were no significant changes in N productivity or NMRT; yet, both slopes were negative mirroring the decline in NUE with increased density (Appendix S2: Table S3; Figure 5e–g). The decline in NUE was concurrent with an observed increase in allocation to foliar N relative to wood N on a per tree basis (Figure 5h). For the most part, when tree-level metrics were calculated on a per unit area, there were few C or N cycling metrics that were significantly related to tree density (Appendix S2: Table S3).

Generate image description

FIGURE 5. Tree-level N cycling traits in relation to increasing tree density. Points are raw data, and the line is the significant linear mixed effects model fit with the 95% confidence interval shaded in gray. A significant linear fit was not observed for the data presented in panels (e) and (f). Model fits for panels (a)–(d) are plotted without raw data in Appendix S3: Figure S3 for ease of viewing. Darkness of green scales from low to high density. Total N is sum of tree-level N in stem, old and new branches, and peak biomass foliage. Total ANPP N is the tree-level sum of N in woody ANPP and foliage. N uptake is the tree-level ANPP N−N lost in September foliage. N resorbed is the tree-level foliar N pool in July−the foliar N pool in September. NMRT is the tree-level mean residence time of N, calculated as the total N pool/N uptake. N productivity is the tree-level total ANPP biomass pool/total ANPP N pool. NUE is the tree-level N use efficiency calculated as the NMRT × N productivity. Foliar:woody ANPP N is the tree-level ratio of N in July foliage to woody ANPP.

Stand-level metrics of C and N cycling almost all showed positive quadratic relationships with tree density with an increase up to ~2.5 trees m⁻² and then a plateau or decline at high densities (Appendix S2: Table S3; Figure 6a–d,h). Of all the stand N cycling variables, only the proportion resorption of N showed no relationship with tree density (Appendix S2: Table S3). Stand-level biomass, productivity, N pools, N uptake, and resorption all increased with tree density up to ~2.5 trees m⁻² and then plateaued or declined, while NMRT, N productivity, and NUE declined with density (Figure 6e–g; Appendix S2: Table S3). The total N pool, ANPP N, N uptake, and N resorption (Figure 6a–d) showed contrasting relationships with tree density at the stand (positive quadratic; Figure 6) and individual tree scales (negative; Figure 5a–d). Similar to observations at the tree level, at the stand level we observed an increase in the ratio of foliage to wood pools as density increased (Figure 6h) that plateaued or declined at ~2.5 trees m⁻² (adjusted R² = 0.47, p < 0.001).

Generate image description

FIGURE 6. Stand-level N cycling characteristics in relation to increasing tree density. Points are raw data, and the line is the linear mixed effects model fit with the 95% confidence interval shaded in gray. Most model fits were improved by employing a quadratic function. Darkness of green scales from low to high density. Total N is the stand-level sum of tree-level N in stem, old and new branch, and peak biomass foliage per unit area. Total ANPP N is the stand-level sum of tree-level N in woody ANPP and foliage per unit area. N uptake is the stand-level sum of tree-level N uptake (ANPP N−N lost in September foliage) per unit area. N resorbed is the stand-level sum of the tree-level foliar N pool in July−the foliar N pool in September per unit area. NMRT is mean residence time of N calculated as the total N pool/N uptake per unit area. N productivity is the stand-level total ANPP biomass pool/total ANPP N pool per unit area. NUE is the stand-level N use efficiency calculated as the NMRT × N productivity. Foliar:woody ANPP N is the mean stand-level ratio of N in July foliage to woody ANPP.

DISCUSSION

We explored how aboveground and belowground C and N cycling metrics were related to increased tree density to improve our understanding of the relationship between N cycling and TTE dynamics with climate warming. Overall, we did not observe consistent evidence of reduced N supply and greater N limitation with increased tree density through the evaluation of soil, tree, and stand aboveground and belowground characteristics. Instead, we observed that as tree density and forest productivity increased, rapid turnover of N due to allocation to short-lived tissues and N uptake supported growth. Increased stand-level allocation to larch roots likely supported the increased N uptake and stand productivity. We did not observe a density effect on resin-sorbed soil N. The contrast between greater stand-level N uptake and no change in resin-sorbed N suggests organic forms of N may be important in supporting productivity. Foliar N traits of individual trees did not indicate reduced N availability with increased density (i.e., no decline in N concentration or δ¹⁵N or increase in resorption). In general, tree N and biomass pools, uptake, and resorption declined on a per tree basis but increased on the stand level up to a density threshold around 2.5 trees m⁻². Opposing processes of decreased NUE and increased N uptake along with increased allocation of N to foliage occurred as stand productivity saturated with increased density. A change in allocation to short-lived tissues with relatively high N content (foliage) reduced the residence time of N in biomass and overall NUE. In contrast to our expectations, foliar traits and allocation patterns suggest competition for aboveground resources may impact C and N cycling more so than N competition.

Does N availability vary with stand density?

Tree density was negatively related to active layer thickness, but we observed no relationship between tree density and SOL thickness or resin-sorbed labile N. Our findings support those observed in stands of a congeneric, L. kaempferi where the experimental manipulation of density through thinning did not result in changes in inorganic N concentrations measured using ion exchange resins (Son et al., 2004). Availability of N in soils reflects both N supply and demand (Craine et al., 2015). Across the density gradient, variation in N availability may be masked by three processes: (1) N leaching in low-density stands, where deeper active layers provide more soil volume for vertical loss to occur; (2) greater N demand and uptake of the available N pool by more individuals in high-density stands, which is supported by our stand-level observations; and (3) interspecific competition in low-density stands where there is high abundance of tall shrubs (Paulson et al., 2021), which are ectomycorrhizal like larch trees, countered by intraspecific competition in high-density stands both of which exert strong impacts on soil N. Furthermore, there are methodological limitations to measuring N availability with resins due to their limited ion sink and measurement solely of inorganic forms. Our observation of greater stand-level N uptake with increased density suggests that organic N may be an important source supporting productivity across our study system as has been suggested for other Siberian larch forests (Tokuchi et al., 2010). Further assessment of variation in N supply and turnover is necessary and could be more comprehensively addressed by evaluating the bulk soil C:N, extractable N, N mineralization rates, and organic N pools.

Do stand-level root characteristics indicate covariation of N availability and density?

Increased uptake requires either increased N availability or increased allocation to N acquisition. We found no evidence for density effects on soil N availability at least for labile forms of N. The fine root biomass of individual trees declined with density along with tree-level size and productivity (i.e., small trees have less roots). However, at the stand level, larch coarse-root density and the proportion of larch fine-root density increased with tree density suggestive of investment in acquiring limited soil resources (Bloom et al., 1985; Nadelhoffer, 2000). Greater root density reduces the likelihood that N in soil solution will be immobilized before reaching the root (McMurtrie et al., 2012) and provides greater foraging surface area within the soil column. Rhizodeposition can stimulate decomposition of soil organic matter and thus increase N available for uptake particularly in N-limited systems (Dijkstra et al., 2013). Furthermore, larch are obligately ectomycorrhizal, and as root density increases allocation to mycorrhizal taxa with exploration types that are particularly adept at nutrient acquisition from labile forms, relatively low C costs, low immobilization in fungal biomass, and high transfer of N to hosts are expected to dominate (Hobbie & Agerer, 2010; Peay et al., 2011). With these mechanisms potentially at play, greater N uptake could support stand-level increases in productivity.

Larch root tissue chemistry is likely reflective of soil N availability and demand. Fine roots in the organic horizon had lower C:N and more depleted δ¹⁵N than fine roots in mineral soil, suggesting these roots may be less reliant on mycorrhizal fungi (Hobbie & Colpaert, 2003) or they are taking up more depleted forms of N (Templer et al., 2007). This contrasts, however, with our observations of greater resin-sorbed NO₃⁻-N, a relatively depleted N source, in mineral but not organic soils. The greater NO₃⁻-N in the mineral soils may reflect lower demand by larch from this soil layer and thus a greater amount in the soil, given that larch and other species were not observed to have appreciable amounts of fine roots in mineral soils and N uptake from mineral soils is lower than that from organic soils (Hewitt, Miller, et al., 2022a).

Does variation in foliar traits reflect N limitation?

Under conditions of high competition and scarce soil nutrients, we expected to see foliar traits of individual trees that indicate greater conservation of nutrients and N limitation as density increased, that is, a decline in N concentration and δ¹⁵N and an increase in resorption. Yet, there was no change in foliar N concentration or δ¹⁵N at peak biomass nor resorption of N. We observed N resorption (~73%) similar to resorption by congeneric L. kaempferi seedlings (76%) under nutrient starvation conditions (Yan et al., 2018) and higher than the global mean (62%; Vergutz et al., 2012). Furthermore, during senescence in September, foliar δ¹⁵N decreased with increasing tree density and greater proportion of N resorbed. A declining trend in δ¹⁵N from pre- to post-abscission has been documented for congeneric L. laricina (Chapin III & Kedrowski, 1983) and L. sibirica needles with a similar magnitude of resorption (Hayashi et al., 2018), and across many other species (Enta et al., 2020). These patterns in foliar δ¹⁵N coupled with the high resorption levels confirm that these TTE forests are low-N environments but indicate that N limitation does not increase with density, potentially because of shifts in allocation and turnover of N that support productivity.

There were strong correlations between foliar δ¹⁵N and δ¹³C in both July and September, suggesting a coupling of larch nutrition and parameters influencing photosynthesis throughout the growing season. However, in contrast to our expectations, it was the foliar C content and δ¹³C but not the N traits that showed the greatest response to increased tree density, thus indicating the importance of non-nutrient resources on foliar stoichiometry. Foliar C content and δ¹³C are related to light and water availability. Across the gradient, light intensity declined but soil moisture did not vary systematically (Loranty & Alexander, 2021). Thus, these foliar C patterns could reflect a greater proportion of C assimilation from sunlit leaves that show less δ¹³C discrimination in low-density stands than shaded leaves in high-density stands (Francey et al., 1985). The needles of trees in high-density stands within our study gradient have higher specific leaf area (Kropp et al., 2019) than those in low-density stands, which is a trait of leaves that are shaded (Liu et al., 2016), further indicating that shading may be an important factor affecting foliar traits as density increases. The range of foliar July δ¹³C signatures we observed (−29.1 to −32.3‰) is more depleted than those reported for eight larch species (−25.4 to −29.5‰; Kloeppel et al., 1998) and declined in July δ¹³C with increased tree density countering conclusions of increased water stress along the density gradient (Kropp et al., 2019; Walker et al., 2021). Despite weak trends between density and most foliar N traits, the relationship between C and N isotope signatures at both peak biomass and senescence implies a coupling of C and N that might be more strongly driven by non-nutrient competition with greater tree density.

Does variation in stand-level C and N cycling characteristics reflect N limitation?

Nitrogen cycling characteristics differed in response to increased tree density by ecological level, from an individual tree to the ecosystem. On a per tree basis, N pools, N uptake, and N resorption along with root biomass were all negatively related to tree density, tracking density-driven declines in individual tree productivity (Walker et al., 2021). At the stand level, greater production associated with higher tree density (Walker et al., 2021) was coupled with greater N storage, uptake, and resorption, but was curtailed around 2.5 trees m⁻². The plateau and/or decline in stand N metrics at the highest densities indicated a limit to ecosystem C and N accumulation, demand, and resorption. The stasis in N cycling characteristics at high densities mirrors the low responsiveness of growth to climate warming at moderate and high densities, despite a positive climate–growth response at low densities (Walker et al., 2021). Thus, density may influence the degree of response of C and N cycling characteristics to future climate warming.

We expected increased productivity to be associated with increased NUE or greater N uptake. Increased NUE occurs due to either increased biomass production per unit N, that is, N productivity, or increased residence time of the N in biomass. Surprisingly, we observed a consistent decrease in NUE on the ecosystem and tree level with increased density concurrent with no change in soil N and at the stand level a reduction in N residence time and N productivity. The shift in NUE and NMRT was attributed to the increased foliage per unit wood at both the individual tree and stand level. Our findings support those from other N-limited systems where increased N allocation to pools with high N content and rapid turnover, that is, fine roots and leaves, resulted in shorter residence time and underscored the importance of N uptake to increased productivity (Finzi et al., 2007).

Implications for TTE advance

Landscape modeling experiments highlight both inconsistencies in the predictions pertaining to the distribution of larch across Siberia and vulnerabilities of larch to future climate warming and associated disturbances (He et al., 2017; Shuman et al., 2011; Tchebakova et al., 2009). Dynamic vegetation models that represent C–N coupling have shown incongruence between modeled and mapped TTE location in eastern Siberia (Wårlind et al., 2014). Our results suggest that increases in stand-level productivity are supported by greater N uptake likely due to allocation to fast cycling tissues, that is, greater foliage: wood and greater root density, and thus high turnover of N. Regionally applied individual-based models (Kruse et al., 2016; Shuman & Shugart, 2009) could be parameterized to account for changes in allocation and test the impacts of density on forecasts of C and N cycling and the distribution of larch forests with future warming across Siberia. Exploration of constraints on productivity including potential light limitation at the highest densities needs to be further investigated to improve forecasting of the effects of infilling on N cycling and forest productivity at TTE under future scenarios of warming.

ACKNOWLEDGMENTS

This research was supported by the Office of Polar Programs Arctic System Science grants to Michelle C. Mack (1545558), Michelle C. Mack and Rebecca E. Hewitt (1708344), and Heather D. Alexander (1304040 and 1708307), Susan M. Natali (1304007), and Michael M. Loranty (1304464 and 1623764). We thank Melissa Boyd and Samantha Miller for help in the field and laboratory.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

DATA AVAILABILITY STATEMENT

Data (Alexander et al., 2021a, 2021b; Hewitt et al., 2018; Hewitt, Alexander, et al., 2022a; Hewitt, Alexander, et al., 2022b; Hewitt, Miller, et al., 2022b) are available from the Arctic Data Center: https://doi.org/10.18739/A2C24QP3W, https://doi.org/10.18739/A2GT5FG86, https://doi.org/10.18739/A2GX44W2M, https://doi.org/10.18739/A2C824G1H, https://doi.org/10.18739/A27H1DN87, and https://doi.org/10.18739/A23R0PV43.