Soil microenvironmental variation drives below‐ground trait variation and interacts with macroclimate to structure above‐ground trait variation of arctic shrubs

Intraspecific trait variation can influence plant performance in different environments and may thereby determine the ability of individual plants to respond to climate change. However, our understanding of its patterns and environmental drivers across different spatial scales is incomplete, especially in understudied regions like the Arctic. To fill this knowledge gap, we examined above‐ground and below‐ground traits from three shrub taxa expanding across the tundra biome and evaluated their relationships with multiple microenvironmental and macroclimatic factors. The traits reflected plant size and structure (plant height, leaf area and root to shoot ratio), leaf economics (specific leaf area, nitrogen content), and root economics and collaboration with mycorrhizal fungi (specific root length, root tissue density, nitrogen content, and ectomycorrhizal colonisation intensity). We also measured leaf and root δ15N and leaf δ13C to characterise nitrogen source and acquisition pathways and plant water stress. Traits were measured in replicated plots (N = 135) varying in soil microclimate, thaw depth and organic layer thickness established across five sites spanning a macroclimate gradient in northern Alaska. This hierarchical design allowed us to disentangle the independent and combined effects of fine‐scale and broad‐scale factors on intraspecific trait variation. We found substantial intraspecific variation at fine spatial scales for most traits and less variation along the macroclimate gradient and between shrub taxa. Consistent with these patterns, microenvironmental factors, mainly soil moisture and thaw depth, interacted with macroclimate, mainly climatic water deficit, to structure size‐structural and leaf trait variation. In contrast, most root traits responded additively to thaw depth and macroclimate. Synthesis. Our results demonstrate that above‐ground and below‐ground tundra shrub traits respond differently to microenvironmental and macroclimatic variation. These differing responses contribute to substantial trait variation at fine spatial scales and may decouple above‐ground and below‐ground trait responses to climate change.


| INTRODUC TI ON
Functional traits can determine plant responses to environmental change, as well as the effects of plants on ecosystem function (Lavorel & Garnier, 2002).While many studies have focused on differences in mean functional trait values among species, there is increasing evidence that plants display substantial intraspecific variation in functional traits (Niu et al., 2020;Siefert et al., 2015), which can influence individual performance, interactions among plants, and community structure and dynamics (Albert et al., 2011;Violle et al., 2012;Westerband et al., 2021).Intraspecific trait variation can also influence the ability of individuals to cope with climate change (Henn et al., 2018;Nicotra et al., 2010), potentially mediating shifts in species ranges and ecosystem functioning due to climate warming (Anderson & Gezon, 2015;He et al., 2018;Valladares et al., 2014).
A predictive understanding of the patterns and drivers of intraspecific trait variation is therefore critical for anticipating biodiversity changes and their ecosystem-level consequences.
Environmental heterogeneity is considered a key factor maintaining intraspecific trait variation.There is evidence that intraspecific trait variation increases with increasing environmental heterogeneity (Stotz et al., 2021), saturating when the observed scale encompasses the full range of environmental conditions experienced by a species (Albert et al., 2011).Yet research also shows that patterns of intraspecific trait variation are idiosyncratic and depend on the spatial and ecological scales of investigation, as well as trait identity (Messier et al., 2010;Siefert et al., 2015;Westerband et al., 2021).For example, while climate and elevation have frequently been shown to drive intraspecific trait variation at regional spatial scales (Anderegg et al., 2021;Midolo et al., 2019;Ostonen et al., 2011;Zadworny et al., 2016), recent studies have documented substantial intraspecific trait variation within populations at local spatial scales (Betway et al., 2021;Burton et al., 2017;Kumordzi et al., 2019;Weemstra et al., 2021).This suggests that microenvironmental factors, like local soil properties and microclimate, may be just as or even more important for driving intraspecific trait variation (Kemppinen & Niittynen, 2022;Siefert et al., 2014).Multiple environmental drivers could also operate together or potentially interact across scales to structure intraspecific trait variation because multiple resources limit plant growth (i.e.energy, water, nutrients) and their supply varies at different scales (Chapin et al., 1987;Gleeson & Tilman, 1992).However, intraspecific trait responses to multiple environmental drivers at different scales have not been widely investigated.
Elucidating the patterns and drivers of intraspecific trait variation is particularly important for understanding and anticipating the dynamics of Arctic tundra plant communities.The Arctic is warming at a rate nearly four times faster than the global average (Rantanen et al., 2022), leading to profound shifts in tundra vegetation (Heijmans et al., 2022;Macander et al., 2022).One of the most widely documented of these shifts is the expansion and accelerated growth of deciduous shrubs (Forbes et al., 2010;Tape et al., 2006), particularly within the genera Alnus, Betula and Salix.Although our knowledge of the factors that influence shrub dynamics is growing (Martin et al., 2017;Mekonnen et al., 2021), predicting the responses of shrubs to climate warming across the Arctic remains challenging.Empirical studies demonstrate that different deciduous shrub species are increasing at different rates and in different areas, which could be related to their distinct ecological strategies (Rinas et al., 2017;Young et al., 2016).Trait-based analyses could provide insights into the functional and strategic variation among species and the mechanisms governing their current and future distributions (García Criado et al., 2023;Pollock et al., 2012).For example, higher variation in leaf economic traits like specific leaf area (SLA) and leaf nitrogen content could signal greater versatility in resource acquisition and use strategy which may translate to greater niche breadth, while the environmental drivers of this variation could indicate the environmental conditions that maximise shrub growth and fitness through environmental matching (Ackerly, 2003;Cardou et al., 2022), potentially enabling improved predictions of shrub dynamics.
Previous studies in tundra ecosystems show that intraspecific trait variation accounts for 17%-74% of the total variation for common size-structural (plant height, leaf size) and leaf economics traits (Jónsdóttir et al., 2023;Kemppinen & Niittynen, 2022;Thomas et al., 2020).For most traits, intraspecific trait variation has been found to be highest at local scales (<10 km 2 ), suggesting microenvironmental heterogeneity partly structures plant trait distributions (Thomas et al., 2020).Recent investigations support this, indicating that intraspecific trait variation of several tundra plant species is strongly related to microclimate, including local temperature, soil moisture content and snow conditions (Kemppinen et al., 2021;Kemppinen & Niittynen, 2022).Studies have also found that changes in local thaw depth can lead to variation in plant chemistry by altering nutrient availability and plant uptake (Blume-Werry et al., 2019;Jasinski et al., 2022).At the same time, broad-scale studies show consistent relationships between intraspecific trait variation and macroclimate across the tundra biome (Bjorkman et al., 2018; Kudo spatial scales and may decouple above-ground and below-ground trait responses to climate change.

K E Y W O R D S
active layer depth, Arctic shrub expansion, intraspecific trait variation, leaf traits, microclimate, root traits, scaling, thaw depth, trait-environment relationships et al., 2001).Given evidence that soil moisture can also mediate macroclimate effects on tundra shrub growth and abundance (Ackerman et al., 2017;Elmendorf et al., 2012;Myers-Smith et al., 2015), it is likely that both macroclimate and microenvironmental factors drive shrub intraspecific trait variation.Yet few studies have investigated the combined influence of these drivers in arctic tundra ecosystems.
Most studies evaluating intraspecific trait variation in the tundra have focused on above-ground traits (Baruah et al., 2017;Betway et al., 2021;Bjorkman et al., 2018;Hudson et al., 2011;Thomas et al., 2020).However, intraspecific variation of root traits may be of greater importance for understanding and predicting vegetation responses to climate change because tundra plants allocate a large fraction of their biomass below-ground (Iversen et al., 2015).Spitzer et al. (2023) found that high levels of root trait plasticity were more influential than species turnover for driving communitylevel trait responses to environmental variation in sub-arctic tundra.Additionally, because roots interact more directly with the soil microenvironment, root traits may be more sensitive than aboveground traits to these factors, as has been shown for tundra and other ecosystems (Spitzer et al., 2023;St. Martin & Mallik, 2021;Weemstra et al., 2022).Divergent responses in above-ground and below-ground traits to environmental variation at different scales would challenge our ability to anticipate how climate change will affect plant species and communities.
We aimed to improve understanding of the patterns of intraspecific trait variation in the Arctic by investigating how the sizestructural, leaf and root traits of expanding shrubs respond to environmental variation at different spatial scales.We asked: (1) How do patterns of trait variation differ across sites, within and among taxa and across plots?(2) What are the primary environmental drivers of trait variation across different spatial scales?Given evidence that tundra plants, particularly deciduous shrubs, possess a range of ecological strategies that are poorly represented by plant functional groups (Saccone et al., 2017;Thomas et al., 2019) and influence responses to climate warming (Rinas et al., 2017;Young et al., 2016), we predicted that trait variation would be greatest among taxa, followed by within sites, in alignment with the expected influence of the soil microenvironment.We further predicted that both soil microenvironmental and macroclimatic factors would structure shrub trait variation either additively or interactively, given previous findings that soil moisture can affect tundra shrub growth responses to climate warming (Ackerman et al., 2017;Elmendorf et al., 2012;Myers-Smith et al., 2015).

| Study area and sampling design
Our study was conducted in northern Alaska at five sites spanning a latitudinal gradient from the southern foothills of the Brooks Range to the Arctic Coastal Plain (Figure 1a; latitude: 67.0° N-69.3°N; longitude: 148.7°W-150.3°W).The sites were spaced 36-265 km apart, with the two southern sites closest together.All sites are in the continuous permafrost zone (Jorgenson et al., 2008) and classified as moist acidic tundra.For the 10-year period from 2006 to 2015, summer air temperature ranged from 11.3 to 13.0°C, winter air temperature ranged from −25.6 to −22.7°C and mean annual precipitation ranged from 233 to 507 mm (Climatic Research Unit, CRU TS 4.0 datasets; Harris et al., 2020) across the sites, with macroclimatic conditions generally becoming more conducive to plant growth from north to south (Table S1).
The vegetation communities sampled at each site are low shrub tundra dominated by the tussock-forming sedge Eriophorum vaginatum and located on flat or gentle slopes (2.5-10°).At each site, we established 27 circular plots (1 m diameter) centred on an individual from one of three expanding deciduous shrub genera (N = 9 plots/genera/site): Alnus (A. viridis ssp.fruticosa), Betula (B. nana) and Salix (S. glauca, S. pulchra, S. arbuscoloides or S. richardsonii), respectively alder, birch and willow (N = 135 plots total).The selected shrub was the dominant shrub in each plot, and E. vaginatum was either the most dominant or next most dominant species in each plot (Table S2).Plots were distributed along a microtopographic gradient from water tracks to upland areas to capture a wide range of soil microenvironmental conditions within each site (Figure 1b).To ensure sampling of taxa experiencing comparable conditions, we spatially clustered plots by taxa along the microtopographic gradient.We did this by first identifying healthy, reproductively mature individuals from each of the three taxa that were near each other (≤5 m apart) and then establishing these plot clusters along the gradient, resulting in clusters with comparable microenvironments (Figure S1).We avoided shrubs showing signs of damage, for example, from herbivory by Lepus americanus (snowshoe hare).Clusters of plots were located 10-550 m apart (median = 37 m) and at least 100 m from roads to minimise road effects on edaphic characteristics (Ackerman & Finlay, 2019).The Salix species sampled in each site are listed in Table S3.Because willow species frequently hybridise (FNA, 1993), we examined congeneric trait variation for willow and intraspecific variation for alder and birch.

| Environmental measurements
To characterise macroclimatic variation, we retrieved downscaled (1 km) monthly air temperature and precipitation data from the Climatic Research Unit (CRU TS 4.0 datasets; Harris et al., 2020).
From these data, we calculated the following variables for each site, averaged over the 10-year period from 2006 to 2015: mean summer air temperature (June-August), mean winter temperature (December-February), mean annual temperature, mean summer precipitation (June-August) and mean annual precipitation.Mean annual growing season length and mean annual snowfall (based on monthly snowfall equivalent data; Littell et al., 2018) at each site were calculated from downscaled (771 m × 771 m) CRU TS 3.1 for the 10-year period from 2000 to 2009 (Harris et al., 2014).See http:// www.snap.uaf.edu for climate data and details.We also retrieved actual evapotranspiration and climatic water deficit data for the 30year period from 1975 to 2005 for each site from 60-m resolution gridded surfaces produced by Morrison et al. (2019).These water balance variables reflect the interactive effects of energy and water on plants across multiple spatial scales (Stephenson, 1990).Actual evapotranspiration is the mean annual evaporative water loss given the prevailing water availability at a site.Climatic water deficit is the mean annual amount of evaporative demand that is not met by available water (irrespective of soil water content) and is related to drought stress (Stephenson, 1998).
To characterise soil microenvironmental variation, we measured thaw depth and soil microclimate (temperature and moisture) between mid and late July 2017 and 2018.Each variable was measured at three separate locations within each plot in each year, and measurements were averaged by plot across years.Thaw depth was measured to the nearest centimetre by inserting a graduated metal rod (1 m in length) into the ground until the frozen surface was reached.
Soil temperature was measured with a hand-held digital thermometer (10 cm depth; accuracy 0.4°C; resolution 0.1°C), and soil moisture was measured with a TDR probe (12 cm depth; Field Scout TDR 100, Spectrum Technologies; accuracy 3% volumetric water content (VWC); resolution 0.1% VWC).The soil temperature and moisture measurements reflected relative differences in soil microclimate, and generally captured warmer versus cooler, or drier versus wetter microclimates, across the course of the study.Because soil organic layer thickness can affect thermal dynamics of tundra soils and their response to climate warming (Baughman et al., 2015), we also collected an intact soil core (moss to permafrost; 7.6-cm diam) from each plot in 2017 to determine the thickness of the organic layer as described in Chen et al. (2020).

| Trait measurements
Between 17 July and 23 July 2017, we measured 12 shrub traits in each plot.Sites were sampled from south to north to ensure that measurements were collected at peak biomass.Three traits reflected plant size and structure (plant height, leaf area, root to shoot ratio), two reflected leaf economics (SLA, N content), and four reflected root economics and collaboration with mycorrhizal fungi greater mycorrhizal dependence as a result of discrimination against the heavier 15 N over the lighter 14 N isotope during N uptake and transfer from mycorrhizal fungi to host plants (Craine et al., 2015).
Leaf δ 13 C provides an integrated measure of plant water stress, with more positive δ 13 C values indicating water stress as a result of stomatal closure and the corresponding reduction in discrimination against the heavier 13 C isotope that occurs during diffusion of CO 2 through the stomata (Farquhar et al., 1989).The procedures for trait measurements followed standard protocols as described in Perez-Harguindeguy et al. ( 2013) and Freschet et al. (2021).
From each target shrub, we collected 10 sun-leaves for leaf trait measurements and excavated at least six intact root branches (coarse + fine roots) from the upper 20 cm of soil as described in Chen et al. (2020).We also recorded the height (cm) of each target shrub (ground surface to highest point of the plant) and basal diameter before root excavation to calculate above-ground biomass (Berner et al., 2015).Leaf and root samples were stored in coolers in zip-lock bags with dampened paper towels before returning to the lab, where they were kept at 4°C until processing.
In the lab, we removed the petioles from fresh leaves, scanned the leaves (400 dpi) using a flatbed scanner and then dried them at 55°C for 72 h.Leaf area was measured from scans (cm 2 ) using ImageJ software (Schneider et al., 2012) and SLA (cm 2 leaf dry mass g −1 ) was computed by dividing the leaf area by oven-dry mass.
Roots were washed with DI water and a subsample of approximately 15 absorptive roots from each plot was selected for additional analysis.We defined the first three root orders as absorptive roots (with the first order being the most distal; sensu McCormack et al., 2015).Roots were held in place between thin acrylic sheets and scanned at 400 dpi; scans were then analysed with WinRHIZO (Regent Instruments, Quebec, Canada) to determine root diameter (RD) and total root length, which in turn was used to calculate SRL (total root length/total root dry mass, m g −1 ) after oven-drying the samples.Root tissue density was calculated from RD and SRL, assuming the root has a cylindrical shape, as: root tissue density = 4/ (pi × RD 2 × SRL, g cm −3 ).We also measured the biomass of live roots (absorptive + non-absorptive roots (woody)) in the soil cores collected from each plot (Chen et al., 2020) and calculated the ratio of total root biomass to above-ground shrub biomass to evaluate allocation patterns.Another subset of fresh shrub root sample was preserved in 70% (v/v) ethanol for later determination of the percentage of root tips that were colonised by ectomycorrhizal fungi by inspection of root tip morphology, colour, architecture and anatomy as described in Chen et al. (2020).
Dried leaf and absorptive root samples were ground into a fine powder and then analysed for total N, δ 15 N and δ 13 C (leaves only) on a Costech 4010 CHNSO Elemental Analyser (Costech Analytical Technologies Inc., Valencia, CA, USA) interfaced with an isotope ratio mass spectrometer (Thermo Fisher Delta V Advantage, Fisher Scientific) at the University of Illinois.Analytical error was <1% for total N and ±0.2‰ for δ 15 N and δ 13 C.All environmental and trait data used in this study are publicly available (Fraterrigo & Chen, 2023a, 2023b).Additionally, all fieldwork performed to collect the vegetation and microenvironmental data was permitted though the Alaska Department of Natural Resources (#LAS31709) and the Bureau of Land Management (#F97344).

| Data analysis
All statistical analyses were performed in the R statistical software version 4.1.3(R Core Team, 2021).To quantify trait variation, we calculated the coefficient of variation of each trait (trait standard deviation/trait mean × 100%) for each taxa and across taxa.We also partitioned total trait variance across nested spatial scales and biological levels using linear mixed effects models and the nlme package (Pinheiro et al., 2023).We calculated the percentage of trait variance explained by trait differences between taxa, within taxa between sites, and within taxa and sites.We specified trait as the dependent variable and taxonomic identity nested in site as a random effect in the model.The residual variance is explained by trait differences between taxonomically identical plants growing in different plots at the same site (Albert et al., 2010).
To determine the environmental drivers of trait variation, we used Bayesian model selection and multimodel inference, which incorporate uncertainty into predictions and improve model performance (Hoeting et al., 1999).Our approach consisted of three steps: (1) fitting a set of candidate hierarchical linear models with different soil microenvironment (soil temperature, soil moisture, thaw depth, and organic layer depth) and macroclimate predictors for each trait, (2) identifying the highest ranked model across all combinations of predictors, and (3) evaluating the relative importance of predictors and calculating model averaged parameter estimates.We briefly describe these steps below and provide a detailed description in Methods S1.
For each individual trait, we first fit a set of candidate Bayesian hierarchical linear models using Stan (Carpenter et al., 2017;Stan Development Team, 2020) called from R through the brms package (Bürkner, 2017).Candidate models included one of the following as fixed effects: a single soil microenvironment variable, a single macroclimate variable, the additive interaction of single soil microenvironment and macroclimate variables, or the additive and multiplicative interaction of single soil microenvironment and macroclimate variables.The modelled microenvironment variables represented conditions at the plot scale, and the macroclimate variables represented conditions at the site scale.After excluding macroclimate variables that were strongly correlated (|r| > 0.95) (Figure S2), we fit 68 candidate models (including a null model) for each individual trait.All models included taxa as a fixed effect, and taxa nested in site as a random effect, to account for the spatial and taxonomic structure of the data.Continuous predictor variables were scaled by the standard deviation and centred on zero to eliminate scaling issues prior to analysis.Height, leaf area and leaf N content were log-transformed to meet assumptions of normality.We assigned all coefficients noninformative priors to minimise the influence of the prior on the posterior distribution (Gelman, 2006) and estimated the models using a Markov chain Monte Carlo algorithm with four chains run for 3000 iterations after a burn-in of 1000 iterations.Model convergence was confirmed by ensuring that all Gelman-Rubin convergence diagnostic ( R) values were less than 1.1 (Gelman & Rubin, 1992).
To identify the highest ranked model for each trait, we used cross-validation (Hooten & Hobbs, 2015).Cross-validation estimates the average predictive performance of a model by repeatedly fitting and testing the model using splits of the available data into training and test (i.e.out-of-sample) sets (Yates et al., 2023).We used leave-one-out cross-validation (Vehtari et al., 2017) to estimate each model's expected log predictive density, a measure of the model's out-of-sample predictive accuracy, where larger values indicate higher accuracy.We also calculated the difference in expected log predictive density between each model and the highest ranked model.Models differing by <1 were considered competing models.
To combine inference from competing models, we used pseudo Bayesian model averaging (pseudo-BMA; Yao et al., 2018).We first calculated pseudo-BMA weights that are proportional to each model's expected log predictive density and adjusted them using Bayesian bootstrapping, which reduces their sensitivity to the prior specification.Next, we calculated the relative importance of each predictor by summing the adjusted pseudo-BMA weights for each competing model that contained the predictor.Similar to relative importance values computed with the Akaike information criterion (Burnham & Anderson, 2002), these values are interpreted as the proportion of evidence that the predictor is included in the best model (Li & Dunson, 2020).Finally, we drew 4000 samples from the pseudo-BMA posterior to calculate an average estimate, SE, and 95% credible intervals for each parameter.Estimates were averaged over models where the parameter appeared, so they are conditional on the parameter value not being zero and indicate effect size.

| RE SULTS
Most of the variance in deciduous shrub traits was found among individuals within the same taxa (Figure 2).Trait variation within species (alder and birch) or genus (willow) accounted for 61%-90% of the variance in height, root: shoot ratio, SLA, leaf N content, leaf δ 15 N, leaf δ 13 C, root tissue density, root N content, root δ 15 N and mycorrhizal colonisation.For all these traits except leaf N content, most variance (45%-83% of total variance) was found within sites ('ITV within sites' in Figure 2).A relatively small proportion (<1%-24%) of the variance was found between sites among individuals from the same taxa ('ITV between sites' in Figure 2).For leaf N content, variance was partitioned nearly equally among the between (43%) and within site (40%) scales, followed by taxonomic identity (17%).

F I G U R E 2
Variance structure of shrub traits across different ecological scales.ITV, intraspecific trait variation for alder (Alnus viridis ssp.fruticosa) and birch (Betula nana) and congeneric trait variation for willow (Salix); Root:shoot, ratio of root to shoot biomass; SLA, specific leaf area; SRL, specific root length.
In contrast, most variance in leaf area was found among taxa (85%), and the variance in SRL was partitioned between taxonomic identity (49%) and within sites among individuals of the same taxa (51%).
Trait coefficients of variation (CV) tended to be higher for sizestructural and root traits than leaf traits (Figure S3).Across taxa, we observed the highest CVs for root:shoot ratio and root δ 15 N.
The traits with the lowest CVs were leaf δ 13 C (mean = 5%) and SLA (mean = 17%).Considering individual taxa, alder had a higher CV than birch and willow for root:shoot ratio, while willow had a higher CV for mycorrhizal colonisation, and root and leaf δ 15 N.
For 11 of the 12 investigated traits, the highest ranked model included both soil microenvironment and macroclimate parameters (Table 1).Root δ 15 N was the only trait for which the highest ranked model included only a soil microenvironment parameter.Averaged across competing models (Table S4), microenvironment parameters had high relative importance values, indicating strong effects on trait variation (Tables S5-S16).Soil moisture had the highest importance value for leaf economic traits, shrub height, leaf area and root tissue density.Thaw depth had the highest importance value for SRL, root δ 15 N, mycorrhizal colonisation and root:shoot ratio.For root N content, organic layer thickness was the parameter with the highest importance value.
Although the three shrub taxa had different trait values, they showed comparable responses to environmental variation, as indicated by similar slopes (Figures S4-S6).Intraspecific patterns in size-structural traits were largely explained by the non-additive interaction of soil microenvironment and macroclimate parameters (Table 1; Figure 3).Sites with higher water deficits had shrubs of taller stature, but height increased less rapidly with soil moisture in these sites compared to sites with lower water deficits (Table S5).In contrast, root:shoot ratio decreased with thaw depth in sites with higher water deficits, indicating a faster increase in above-ground relative to below-ground biomass as the active layer thickens (Table S7).Leaf area increased with soil temperature in sites with warmer summers, especially alder and willow shrubs, but decreased with soil temperature in sites with cooler summers (Figure S5; Table S6).
Intraspecific patterns in leaf traits were also largely explained by the non-additive interaction of soil microenvironment and macroclimate parameters, specifically soil moisture and climatic water deficit (Table 1; Figure 3).In sites with higher water deficits, leaf economic traits (leaf N content and SLA) increased rapidly with soil moisture, consistent with a shift toward a more resource-acquisitive strategy (Tables S8 and S10).In contrast, leaf N content and SLA decreased with soil moisture in sites with lower water deficits.Leaf 15 N showed a similar pattern, increasing with thaw depth in sites with higher water deficits but not lower water deficits; however, these parameters had credible intervals that crossed zero, indicating high uncertainty (Table S11).Although winter temperature had a lower relative importance value across competing models, it was positively related to leaf 15 N, and its credible interval did not cross zero.An indicator of plant water stress, leaf δ 13 C was the only above-ground trait for which we did not find evidence of an interaction between soil microenvironment and macroclimate parameters (Table S9).Leaf δ 13 C was consistently higher in sites with higher mean annual temperatures and decreased with soil moisture at a similar rate across all sites.
While root traits responded to variation in the soil microenvironment and macroclimate (Figure 4), effects were mostly additive (Table 1; Table S4), with the average relative importance of nonadditive (multiplicative) interactions substantially lower than single soil microenvironment parameters across all root traits (0.05 vs. 0.43, resp.; Tables S12-S16).The exception was SRL, which was interactively affected by summer temperature and thaw depth (Table S12).SRL was higher in sites with warmer summers and decreased more rapidly with thaw depth in sites with cooler summers (Figure S6).By contrast, mycorrhizal colonisation and root δ 15 N increased with thaw depth at a similar rate across all sites (Tables S15 and S16).For

TA B L E 1
The highest ranked model for each trait investigated.All models include taxa as a fixed effect and taxa nested in site as random effect to account for the spatial and taxonomic structure of the data.Bolded terms indicate predictors for which the 95% credible interval did not cross zero based on pseudo Bayesian model averaging over models where the parameter appears (i.e.conditional averaging).R 2 c is the conditional R 2 , and R 2 m is the marginal R 2 .The intraclass correlation coefficient (ICC) indicates the level of similarity among shrubs with the same taxonomic identity within a site.See Table S4 for the full list of competing models, and Tables S5-S16 for relative importance values, model-averaged parameter estimates, estimated error and credible intervals.root N content, organic layer depth was the highest ranked parameter but had a credible interval that crossed zero, whereas root N decreased with mean annual and mean summer precipitation, snowfall, and mean annual temperature (Table S14).For root tissue density, all parameters had credible intervals that crossed zero, indicating high uncertainty (Table S13).

| Overview
We found substantial trait variation among shrub individuals of the same taxa within sites and less variation along the latitudinal gradient and between shrub taxa.Consistent with these patterns, soil microenvironment and macroclimate interactively affected sizestructural and leaf trait patterns.Root traits, in contrast, were additively affected by these factors.Our results demonstrate that arctic shrub trait patterns are heterogeneous at fine spatial scales and strongly structured by a combination of microenvironmental and macroclimatic factors but in different ways.These differential effects of soil microenvironment and macroclimate further suggest that above-ground and below-ground traits will show divergent responses to climate change.

| Magnitude of trait variation at different scales
Root traits generally varied more than size-structural and leaf traits, and most of the variance was found within sites among individuals from the same taxa.This result aligns with other studies that have found substantial intraspecific variation in root traits at local scales (Kumordzi et al., 2019;Spitzer et al., 2023;Weemstra et al., 2021) and suggests tundra shrub roots respond strongly to soil microenvironmental heterogeneity (Hodge, 2004).Interestingly, willow had the highest CV for mycorrhizal colonisation, as well as root and F I G U R E 3 Variation in size and leaf traits across all shrub taxa in relation to the primary soil microenvironment and macroclimate variables identified.Macroclimate variables were modelled as continuous but are shown as two discrete classes, separated by the median value, to improve visualisation.Sites with high water deficit have mean annual water deficits >7.44 mm, while sites with low water deficit have mean annual water deficits ≤7.44 mm.Sites that are warm annually have mean annual air temperatures >−7.15°C, while sites that are cool annually have mean annual air temperatures ≤−7.15°C.Points represent predictions from the best fitting model.In all panels, a single asterisk indicates that the 95% credible interval on the slope of the microenvironment-trait relationship did not overlap zero.Two asterisks indicate that the macroclimate × microenvironment interaction term did not overlap zero.Transparent ribbons are 95% credible intervals for model mean predictions.Corresponding model statistics are presented in Table 1 and Tables S5-S11.Trait-environment relationships by taxa are presented in Figure S5.Note the log scale for height and leaf nitrogen content.leaf δ 15 N.Because δ 15 N reflects dependence on mycorrhizal fungi for N (Craine et al., 2015), these high CVs may indicate contrasting strategies for soil resource acquisition among the willow species we sampled, which could enable willows to occupy a wider range of environments and expand more rapidly across a warming tundra.In contrast, leaf traits like SLA are known to vary more strongly with light and temperature than with soil resource availability (Poorter et al., 2009).In the tundra, deciduous shrubs are unlikely to be light limited because of their tall stature compared to other tundra growth forms, which may be why relatively less leaf trait variance was found within sites.

| Determinants of above-and below-ground trait variation
Previous studies indicate that climate and soil properties at the macro scale jointly influence global and regional patterns of sizestructural and leaf trait variation but also have limited ability to explain community-level trait patterns (Bruelheide et al., 2018;Joswig et al., 2022;Maire et al., 2015;Ordonez et al., 2009).These findings suggest a strong role for local environmental conditions in structuring trait variation.Using a multiscale approach, we found that interactions between soil microenvironmental factors, mainly soil moisture and thaw depth, and climatic water deficit were important predictors of intraspecific variation in shrub size-structural and leaf economic traits in the Alaskan tundra.A synthesis of studies from across the tundra biome showed that community-weighted mean size and leaf trait responses to summer warming varied with sitelevel moisture availability, with leaf area and leaf economic traits leaf N and SLA increasing more rapidly with warming in wet compared to dry sites (Bjorkman et al., 2018).Our results expand on these findings by demonstrating that soil moisture and other soil microenvironmental factors like thaw depth interact with macroclimate factors at finer scales than previously documented and influence species-and genus-level trait responses in addition to communitylevel responses.The strong climate-contingent responses of leaf economic traits to increased soil moisture observed in our study also underscore the importance of phenotypic shifts for maximising resource acquisition in response to both changing microsite and regional conditions.
Similar to other studies, we found taller shrubs in warmer sites with higher water deficits (Bjorkman et al., 2018;Kemppinen & Niittynen, 2022;Walker, 1987).Increased plant size with temperature can be attributed to colder temperatures restricting growth rates (Hudson et al., 2011).However, our study revealed that shrub height responses to climatic water deficit varied within sites, increasing more gradually with soil moisture in warmer sites compared to cooler sites with lower water deficits.High climatic water deficit corresponds to high atmospheric demand for water (Stephenson, 1990), which can limit plant growth by reducing stomatal conductance (Oren et al., 1999;Restaino et al., 2016).Although increased soil moisture can buffer plant water stress associated with atmospheric dryness (Kropp et al., 2017), recent studies indicate atmospheric demand decreases surface conductance and evapotranspiration to a greater extent than soil moisture (Fu et al., 2022;Novick et al., 2016).
The weaker response of shrub height to soil moisture in warmer sites may therefore reflect the overriding effects of higher atmospheric water demand.The leaf 13 C patterns we observed support this hypothesis.Leaf 13 C values declined with soil moisture but were higher in warmer sites, suggesting that plants were consistently experiencing greater water stress in these areas.
In contrast, SLA and leaf N of the different shrub taxa increased with soil moisture in sites with higher water deficits, suggesting that local moisture availability constrains the ability of arctic shrubs to adjust to warming.An increase in SLA signifies lower investment in leaf construction costs to maximise leaf surface area and light capture, and is associated with a shift towards a resource acquisitive strategy (Poorter et al., 2009;Wright et al., 2004).It is therefore expected to increase with temperature.However, leaf economic traits have been found to both increase and decrease along latitudinal and elevation gradients in arctic tundra ecosystems (Betway et al., 2021;Kudo et al., 2001).For example, Betway et al. (2021) found that SLA of Salix pulchra increased from north to south, while Betula nana decreased.One potential explanation for these discrepancies is that the communities or sites where the species occur differ in soil moisture (Baruah et al., 2017;Betway et al., 2021).Overall, our results highlight the need to account for microenvironmental conditions such as moisture and snow depth and cover that can limit plant response to climate change (Kemppinen & Niittynen, 2022).
While root traits responded to both microenvironmental and macroclimatic factors, effects tended to be additive and dominated by microenvironmental factors.Thaw depth strongly predicted SRL, mycorrhizal colonisation, and root 15 N, traits that are linked to resource acquisition strategy (Bergmann et al., 2020;Chen et al., 2020;Weigelt et al., 2021).The contrasting relationships of thaw depth with SRL (negative) and mycorrhizal colonisation (positive) suggest greater dependence of shrubs on mycorrhizal fungi for nutrient acquisition with increasing thaw depth.However, root 15 N also increased with thaw depth, indicating root uptake of inorganic N increased simultaneously (Craine et al., 2015).Other studies have found inorganic N availability and plant N uptake increase with thaw depth in tundra ecosystems, the latter pattern generally attributed to increasing root length and rooting depth (Blume-Werry et al., 2019;Hewitt et al., 2019;Keuper et al., 2017;Pedersen et al., 2020).Thus while mycorrhizal fungi may have contributed to N scavenging from deep soils (Hewitt et al., 2020), deeper rooting by shrubs may also account for the root 15 N patterns observed in our study.
Our results support the emerging view that fine-scale variation in environmental conditions strongly shape trait distributions (Baruah et al., 2017;Betway et al., 2021;Kemppinen & Niittynen, 2022;Spitzer et al., 2023).Furthermore, our results indicate that climate change will differentially affect above-ground and below-ground trait distributions due to the contrasting influence of the local growing environment on trait responses, which could decouple above-and below-ground traits (Weemstra et al., 2022).
Interestingly, microenvironmental conditions, including soil moisture and temperature, were found to consistently predict community-weighted mean size and leaf traits across a large climate gradient when fixed values for species traits (i.e.species-level trait means) were used (Kemppinen et al., 2021).This suggests that non-additive interactions between factors across spatial scales may more commonly structure intraspecific trait variation than community trait variation.Additional work focusing on intraspecific trait variation is needed to better understand how environmental factors interact across spatial scales to shape vegetation responses.
Overall, the patterns of trait variation we observed in this study broadly align with results from previous research which found substantial spatial heterogeneity in shrub growth responses to warming (Myers-Smith et al., 2015, 2020;Reichle et al., 2018).
Moisture availability has been shown to influence the sensitivity of shrub growth to summer warming, with higher sensitivity where moisture availability is greater (Ackerman et al., 2017;Buchwal et al., 2020;Campbell et al., 2021;Niittynen et al., 2020;von Oppen et al., 2021).In addition to providing greater access to water, areas with higher moisture availability may also support higher rates of nutrient mineralisation (Binkley et al., 1994;Chu & Grogan, 2010), allowing shrubs to better exploit increased air temperatures compared to drier areas.This hypothesis is supported by our leaf N results, which show N concentration increased with increasing moisture in warmer but not cooler sites.In contrast, higher moisture in cooler sites has been shown to inhibit growth, possibly as a result of anoxia during summer or increased cold penetration in the winter (Tape et al., 2012).

F I G U R E 4
Variation in root traits across all shrub taxa in relation to the primary soil microenvironment and macroclimate variables identified.Macroclimate variables were modelled as continuous but are shown as two discrete classes, separated by the median value, to improve visualisation.Sites with warm winters have mean winter air temperatures >−23.9°C, while sites with cold winters have mean winter air temperatures ≤−23.9°C.Sites with warm summers have mean summer air temperatures >11.4°C, while sites with cool summers have mean summer air temperatures ≤11.4°C.Sites that are wet annually have mean annual precipitation >282 mm, while sites that are dry annually have mean annual precipitation ≤282 mm.Sites with wet summers have mean summer precipitation >137 mm, while sites with dry summers have mean summer precipitation ≤137 mm.Points represent predictions from the best fitting model.In all panels, a single asterisk indicates that the 95% credible interval on the slope of the microenvironment-trait relationship did not overlap zero.Two asterisks indicate that the macroclimate × microenvironment interaction term did not overlap zero.Transparent ribbons are 95% credible intervals for model mean predictions.Corresponding model statistics are presented in Table 1 and Tables S12-S16.Trait-environment relationships by taxa are presented in Figure S6.

| Implications for modelling vegetation dynamics
Although the three deciduous shrub taxa we investigated varied widely in their trait values at similar positions along environmental gradients, they generally showed similar responses to microenvironmental and macroclimatic variation, consistent with the plant functional type concept.Many Earth system and ecosystem models that simulate vegetation dynamics use the plant functional type concept to simplify vegetation complexity and its effects on ecosystem processes (van Bodegom et al., 2012;Wullschleger et al., 2014).The functional attributes (traits) of a plant functional type are represented with constant parameter values, which implicitly neglects within-and betweenspecies variation and how they relate to environmental conditions.A newer 'next generation' of dynamic vegetation models focuses on including trait variability (Berzaghi et al., 2020;van Bodegom et al., 2014) and aligns well with our results that indicate variation in trait values among shrub taxa across multiple environmental gradients.Indeed, some recent work has focused on applying this trait-based approach to tundra vegetation, finding a strong relationship between field data and trait-based model outputs compared to more traditional modelling approaches (Zhu et al., 2016).Another recent model application has found that when applying a dynamic vegetation model across the same shrub tundra transect described in this present study, both model sensitivity and uncertainty varied spatially along the transect (Euskirchen et al., 2022).Furthermore, this same study found that different shrub taxa (e.g.dwarf birch shrubs compared to other deciduous shrubs such as willow) were sensitive to different model parameters, including parameters related to leaf area, fine roots, and photosynthesis, thereby highlighting the complexity in accurately forecasting vegetation change in the highly heterogeneous arctic tundra plant communities.

| CON CLUS IONS
Our study reveals large variations in above-and below-ground shrub traits at different spatial scales.Non-additive interactions between soil microenvironment variables and climatic water deficit suggest that above-ground shrub traits will exhibit responses to climate change that are strongly contingent on subsurface environmental conditions, potentially amplifying existing and projected patterns of vegetation heterogeneity at fine spatial scales (Epstein et al., 2004;Lantz et al., 2010).These responses are likely to differ from those of root traits, which are additively affected by these factors.Ultimately, our findings indicate that a better understanding of the soil microenvironment will be needed to accurately predict how tundra plant communities and functional traits will respond to climate change in the Arctic.

F
Locations of the five sampling sites along a latitudinal gradient in northern Alaska (a).Background colour represents mean annual temperature (MAT) for the 10-year period from 2006 to 2015 (Source: CRU TS 4.0 datasets).The upper inset shows the region of Alaska where the sampling transect is located (black line).Soil microenvironmental variables measured in situ (b).Data are presented as violin plots overlaid with box plots.In the box plots, the notches and hinges represent the 25th, 50th (median) and 75th percentiles, the whiskers represent 95% confidence intervals, with outliers shown as points.In the violin plots, the thickness of the polygon corresponds to the local density of the observations.VWC, volumetric water content.(specific root length [SRL], root tissue density, N content, and ectomycorrhizal colonisation intensity).Stable N and C isotope signals in plant tissues integrate plant-environment interactions over long time periods.We therefore measured leaf and root δ 15 N and leaf δ 13 C in each plot.Leaf and root δ 15 N indicate N source and N acquisition pathways, with more negative δ 15 N values indicating