Highlights

  • Impacts of dwarf shrubs on tundra soils were modelled with structural equation models.

  • We investigated summertime soil microclimate and soil organic carbon stocks.

  • Dwarf shrub dominance decreases soil moisture, temperature and organic carbon stocks.

Introduction

Climate change has led to a biome-wide increase in woody plant dominance in the tundra, that is, shrub expansion or shrubification (Myers-Smith and others 2011). Evidence and consequences of the woody plant expansion have been explored with remote sensing (Lantz and others 2013), dendrochronology (Ackerman and others 2017), isotope (Cahoon and others 2016), and experimental studies (Bjorkman and others 2020). Shrub expansion results from favourable environmental conditions, mostly from increasingly warmer and longer growing seasons (Wilson and Nilsson 2009; Hallinger and others 2010; Myers-Smith and Hik 2018). Contrarily, woody plants may feedback to the global climate system through changes in, for example, albedo, but also through changes in soils, such as soil microclimate and carbon stocks (Myers-Smith and others 2011). After all, woody plants are the largest and very common plant life form in the tundra, especially in the sub- and low-Arctic ecosystems, and therefore, extensive changes in their abundance and distribution are likely to lead to cascading effects across the tundra biome and beyond (Swann and others 2010; Cahoon and others 2012).

Overall, plants regulate soil microclimate and the carbon cycle, through for instance, transpiration (Bonfils and others 2012), shading (Myers-Smith and Hik 2013; Loranty and others 2018; Robinson and others 2019), and primary production (McLaren and others 2017) (Figure 1). In the tundra, woody plants intercept rainfall, which, on a large scale, may decrease the overall water input into the soil (Zwieback and others 2019). However, vegetation can also cause opposing effects by shading, which in turn decreases summertime soil temperature (Blok and others 2010; Lantz and others 2013) and evaporation from the soil surface (Bonfils and others 2012). Woody plant expansion has been shown to increase the aboveground carbon storage and cause changes in litter quality and increase biomass production (DeMarco and others 2014; Parker and others 2015), which may further alter the rate of carbon cycling in the tundra (Myers-Smith and Hik 2013). This is globally important as the slow decomposition rate in the cool northern permafrost region has enabled the soils to store a large amount of carbon, roughly 50% of the global belowground organic carbon pool (Hugelius and others 2014).

Figure 1
figure 1

The mechanisms through which woody plants may affect tundra soils. We build a theory based hierarchical model, in which we control for the influence of topography and snow on both woody plants and local soil properties. We expect that the dominance and height of woody plants have a mediating effect on multiple soil properties, namely growing-season soil microclimate and soil organic carbon stocks.

Full size image

Investigations on the ecosystem effects of woody plant expansion are often based on relatively limited data sets. Consequently, the effects of woody plants on soil properties can be challenging to partition from plants themselves or from other environmental factors, such as topography (Crofts and others 2018). In addition, the impacts of shrub expansion have primarily been studied with tall, deciduous shrub species, whilst smaller, albeit more abundant, evergreen dwarf shrubs are understudied, although they differ in litter quality and quantity, and possibly in their impacts on soil microclimate (Vowles and Björk 2019). Nonetheless, evergreen and deciduous dwarf shrubs are responding to climate change (Wilson and Nilsson 2009; Hallinger and others 2010; Maliniemi and others 2018).

Here, we quantify the influence of dwarf shrubs (< 50 cm) on soil and ask “Do woody plants have mediating effects on multiple soil properties in the tundra?”. With “mediating effects” we mean that the woody plants act as moderators on how the macroclimatic and topographic effects on soil conditions are manifested at the local scale, demonstrating the hierarchical nature of these systems. We use 171 plots with in situ measurements and high-resolution remotely sensed data collected systematically across a mountain tundra landscape. We build a structural equation model to test our hypothesis: we expect that despite the strong influence of topography, snow, and the overall vegetation cover, also the low-statured woody plants have separable effects on growing-season soil moisture, soil temperature, and soil organic carbon stocks (Figure 1). We test this with the overall dwarf shrub community and separately for evergreen and deciduous communities. Although woody plants have multiple and potentially contrasting effects on soils, we expect to find subtle, yet direct effects on soil microclimate and carbon stocks (Figure 1).

Materials and Methods

Study Area

We conducted our research in the sub-Arctic mountain tundra of north-western Finland (69° 03′ N 20° 51′ E). On average, the mean annual temperature is − 3.1 °C and annual precipitation 518 mm (1991–2018), first measured at the nearby Saana meteorological station (69.04 N, 20.85 E; 1002 m a.s.l.) and the latter at Kyläkeskus meteorological station (69.04 N, 20.80 E; 480 m a.s.l.), which are located about 1.5 km and 1.0 km from the study site (Finnish Meteorological Institute 2019a, b). The landscape in our study area is characterised by varying topography (Figure 2) and predominantly thin soils, 13.0 cm on average measured using a metal probe (for more details, see Kemppinen and others 2018).

Figure 2
figure 2

Study design in the heterogeneous mountain tundra system. The maps present the fine-scale variation and large gradients in elevation, radiation, topographic wetness index, and topographic position index values across the landscape. White dots represent the 171 study plots (1 m2), from which the in situ measurements were collected. The maps are based on the openly available light detection and ranging-based digital terrain model (2 m × 2 m resolution) provided by the National Land Survey of Finland (2019).

Full size image

In Fennoscandia, the abundance and coverage of dwarf shrubs have increased as a consequence of increasingly warmer summer temperatures and thicker snow cover during winter (Wilson and Nilsson 2009; Hallinger and others 2010; Maliniemi and others 2018). Especially the evergreen Empetrum nigrum has increased significantly during the past decades (Wilson and Nilsson 2009; Maliniemi and others 2018). In our study area, the main vegetation type is dwarf shrub tundra, which forms c. 25–30% of the overall Arctic vegetation (Walker and others 2018; Raynolds and others 2019). The dominant woody plant species are evergreen Vaccinium vitis-idaea and E. nigrum, and deciduous Betula nana and Vaccinium myrtillus.

Study Design

In the tundra, transition zones between habitats create variability in species composition and ecosystem functions (Billings 1973; Fletcher and others 2012). Thus, we selected the location of the 171 study plots (1 m2) using a systematic grid approach (3 km2), which covers a large spatial extent maximising the diversity of habitats as well as the transition zones between them. We excluded plots situated in river channels (with the exception of seasonal meltwater channels) or boulder fields or those exposed to anthropogenic disturbance (that is, trails). In addition to the systematic grid with a 50-m minimum distance between plots, we intentionally situated 25 plots in snow-bank environments and along windswept ridge tops to maximise the snow accumulation gradient in the data. We recorded the locations of the plots using a hand-held Global Navigation Satellite System receiver with an accuracy of up to ≤ 6 cm under optimal conditions (GeoExplorer GeoXH 6000 Series; Trimble Inc., Sunnyvale, CA, USA).

Soil Data

Soil Moisture

We collected the soil moisture data during the 2017 growing season (following Kemppinen and others 2018). We used a hand-held time-domain reflectometry sensor (FieldScout TDR 300; Spectrum Technologies Inc., USA), with a resolution of 0.1% with a reported accuracy of ± 3.0 volumetric water content (VWC%). Moisture was measured from the topsoil (0.0 to 7.5 cm), as the soils in the area are relatively shallow, and therefore, longer measuring probes were not an option.

In our study area, the overall spatial pattern of soil moisture remains similar throughout the growing season (Kemppinen and others 2018), meaning that for instance, depressions are relatively more moist compared to ridges (Billings 1973). However, we repeated the measurements on five occasions (hereafter, campaigns). During each campaign, we took the mean over three measurements per plot to account for possible fine-scale moisture variation within a plot. In subsequent analyses, we used the mean of these five campaigns to represent the spatial pattern of the growing season soil moisture (hereafter, soil moisture).

Soil Temperature

We collected the soil temperature data simultaneously with the soil moisture measurements in 2017. We used a hand-held digital temperature device (VWR-TD11; VWR International, USA), with the resolution of 0.1 °C with a reported accuracy of ± 0.8 °C. Temperature was measured from the topsoil (6.0 to 7.5 cm depth).

We measured temperature once from the centre of each plot and repeated this on the five campaigns. In our analyses, we used the mean over the five campaigns to represent the spatial pattern of growing season soil temperature (hereafter, soil temperature). We corrected the possible effects of the timing of the measurements using data from miniature temperature loggers, which were chiefly installed in the same plots (for details, see Supplementary material 1).

Soil Organic Carbon Stock

We collected soil samples (c. 1 dl) from the close proximity of the plots in August 2016 and 2017. We collected the samples from the organic and mineral soil layers using metal soil core cylinders (4–6 cm Ø, 5–7 cm in height). Organic samples were collected from all plots and mineral samples from a subset of plots. We measured the depth of the organic and mineral soil layers in all plots up to a depth of 80 cm using a metal probe. We took the mean over three depth measurements per layer and per plot to account for possible fine-scale soil depth variation within a plot. Some of our plots were located in tundra meadows (< 10%), in which the topsoil (c. 0–30 cm) can represent a mixture of both organic and mineral layers. Thus, distinguishing the organic layer from the mineral may be reflected in our estimations of the stocks.

First, we analysed the samples at the Laboratory of Geosciences and Geography and the Laboratory of Forest Sciences, University of Helsinki. We used a freeze-drying method to remove soil moisture from the samples. Then, we sieved the mineral soil layer samples through a 2-mm plastic sieve and homogenised the organic soil layer samples by hammering the material into smaller particles.

Secondly, we analysed the bulk density from all organic samples, as well as the total carbon content (C%) or the soil organic matter content (SOM%). Chiefly, we used C% for the majority of the organic layer stock calculations, but we used SOM% for plots that were missing C% data (n = 52). We estimated the bulk density (kg/m3) by dividing the dry weight by the sample volume. We analysed C% using elemental analysers (vario MICRO cube and vario MAX cube; Elementar Analysensysteme GmbH, Germany). We measured SOM% using the loss-on-ignition method according to SFS 3008 (1990). We repeated 17 organic samples in 2016 and 2017, and from these samples we took an average stock across the two years.

In addition, we analysed bulk density from 70 mineral samples and C% from 64 samples, respectively. We used the median bulk density and C% of these samples (830 kg/m3 for bulk density, range 470–1400 kg/m3; 3% for C%, range 0.4–6.5%) to estimate the carbon stock of the mineral layer. We used the median values, as we did not have data on the mineral soil conditions from all plots. We acknowledge the possible uncertainty in our soil organic carbon stock estimates due to the averaging of the mineral samples across plots. However, we argue that in this landscape, the depth of the organic layer and its carbon content are a more critical part of the total stock estimation, as the organic layer stores the majority of the soil carbon compared to the mineral layer.

Thirdly, we calculated the soil organic carbon stocks for the organic layer and mineral layer following Eq. 1 (following Parker and others 2015):

$$ \begin{aligned} {\text{Organic}}\;{\text{layer}}\;{\text{carbon}}\;{\text{stock}} & = {\text{organic}}\;{\text{layer}}\;{\text{C}}\% /{1}00*{\text{organic}}\;{\text{layer}}\;{\text{bulk}}\;{\text{density}}\;\left( {{\text{kg}}/{\text{m}}^{3}} \right)*{\text{organic}}\;{\text{layer}}\;{\text{depth}}\left( {\text{m}} \right),\;{\text{and}} {\text{Mineral}}\;{\text{layer}}\;{\text{carbon}}\;{\text{stock}} & = {\text{mineral}}\;{\text{layer C}}\% /{1}00*{\text{mineral}}\;{\text{layer}}\;{\text{bulk}}\;{\text{density}}\;\left( {{\text{kg}}/{\text{m}}^{3}} \right)*{\text{mineral}}\;{\text{layer}}\;{\text{depth}}\;\left( {\text{m}} \right), \\ \end{aligned} $$
(1)

For plots that were missing organic layer C%, SOM% was used instead of C% in Eq. 1 to calculate the soil organic matter stocks, which we converted into organic layer carbon stocks based on a relationship between C% and SOM% (carbon fraction in the soil organic matter). For the conversion, we used data from plots that had both C% and SOM% data based on the organic layer samples (n =  118). We calculated the carbon fraction in the soil organic matter as 0.54 (R2 = 0.97), similar to Parker and others 2015. We arrived at this fraction by regressing the carbon stock in the organic soil layer using the soil organic matter stock in the organic soil layer without the intercept.

Finally, we summed the organic and mineral layer carbon stocks into the total soil organic carbon stock up to 80 cm (hereafter, soil organic carbon stocks) using Eq. 2:

$$ {\text{Soil}}\;{\text{organic}}\;{\text{carbon}}\;{\text{stock}} = {\text{organic}}\;{\text{layer}}\;{\text{carbon}}\;{\text{stock}} + {\text{mineral}}\;{\text{layer}}\;{\text{carbon}}\;{\text{stock}}, $$
(2)

Vegetation Data

A vascular plant species expert collected the vegetation data from 16 July to 8 August 2017. They estimated the coverage of vascular plant species in each plot. In the analysis, we use the plot-specific absolute coverage of vascular plant species (hereafter, plant coverage) and the relative coverage of woody plant species (hereafter, woody plant dominance), also separately for the evergreen and deciduous woody plant species (hereafter, evergreen dominance and deciduous dominance). Woody plant height was measured as the median height across three measurements of each woody plant species in a plot (hereafter, woody plant height). We also use height data separately for the evergreen and deciduous woody plant species (hereafter, evergreen height and deciduous height). We collected the height data from 23rd to 27th July 2016 (146 plots) and completed it between 16 to 19th July 2017 (25 plots). The study plots were under the same grazing pressure, and the main herbivores in the region are Rangifer tarandus tarandus and Cricetidae sp.

In total, we found 19 woody plant species (for species list, see Figure 3). In the height measurements, we excluded possible inflorescence to avoid any bias, particularly regarding the most ground-dwelling species such as Harrimanella hypnoides and Linnaea borealis. Besides the prostrate Salix species (namely, Salix herbacea and Salix reticulata), we grouped all erect Salix species into one category (Salix spp.) given the large amount of Salix hybridisation. We followed the taxonomy of The Plant List 2013.

Figure 3
figure 3

Woody plants of the dwarf shrub dominated tundra. Woody plant dominance and height data in relation to the overall plant coverage and each other from study plots (n = 171). The coverage in relation to the overall plant coverage is shown as bars (A). The relative coverage grouped by deciduous and evergreen species is presented as a frequency plot (B). The prevalence of the species ranks them from the most common (Vaccinium vitis-idaea) to the rarest (Betula pubescens) woody plant species occurring in the data (C). Species-specific coverage and median height data are shown as box plots (D). In the box plots, the notches and hinges represent the 25th, 50th, and 75th percentiles. In addition, the whiskers represent the 95% percentile intervals and the points represent the outliers. The coverage and height in relation to each other are shown as a density plot (E). The hue intensity indicates the local point density.

Full size image

Snow Data

We measured the snow depth to represent the winter conditions in the study system (Aalto and others 2018). We collected snow depth data from 13 to 17th April 2017, during the season of the maximum snow depth. We measured snow depth once per plot from the centre of the plot using an aluminium probe.

Topographic Data

We obtained the light detection and ranging (LiDAR)-based digital terrain model (DTM) at a 2 m resolution from the open file service of the National Land Survey of Finland (2019). Based on the DTM, we calculated elevation, radiation, topographic wetness index (TWI), and topographic position index (TPI) (Figure 2) (following Kemppinen and others 2018).

Radiation

We calculated the potential incoming solar radiation (kWh/m2) for the three months of the growing season (June, July, and August) using a “Potential incoming solar radiation” tool in SAGA GIS v. 2.3.2. We considered the possible shadow effect of obstructing topographic features using the sky view–factor option (Böhner and Antonić 2009). We calculated the position of the Sun for every fifth day using a 4-h interval. We used the lumped atmosphere option to calculate the atmospheric transmittance.

Topographic Wetness Index

We calculated the topographic wetness index (TWI), which is a proxy for the water flow and accumulation using Eq. 3, where SCA refers to the specific catchment area:

$$ {\text{TWI}} = {\text{ln}}\left( {{\text{SCA}}/{\text{local}}\;{\text{slope}}} \right), $$
(3)

TWI can be used to model the variation in the soil moisture at fine-spatial scales (Böhner and Antonić 2009; Kemppinen and others 2018). The total catchment area (TCA) was calculated from a filled DTM using the multiple flow–direction algorithm (Freeman 1991; Wang and Liu 2006). The SCA was calculated assuming that the flow width equals the grid resolution (2 m) (that is, SCA = TCA/2). The local slope (in radians) was calculated using the algorithm by Zevenbergen and Thorne 1987.

Topographic Position Index

We calculated the topographic position index (TPI) based on the elevation difference (m) between a plot and the surrounding elevation along a given radius (Ågren and others 2014). This describes the position of the plot on a topographic gradient—that is, a plot located on a ridge top (positive values), in a depression (negative), or on a slope or flat ground (close to zero). We calculated TPI using 30 m radii based on an unfilled DTM (following Kemppinen and others 2018).

Statistical Analyses

We use structural equation modelling (SEM). SEM produces a single causal network, in which several response variables and predictors are combined by probabilistic models (Grace and others 2010). This enables the testing of hypothesised causal relationships. A variable can be both a predictor and a response variable, that is, a mediator. Thus, SEM allows the simultaneous evaluation of several causal structures, which are expressed as standardised regression coefficients allowing for the direct comparison of the estimated effects (Lefcheck 2016). We use the piecewiseSEM package in R version 3.5.1 (R Development Core Team 2016; Lefcheck 2016).

First, we log-transform the plant coverage, woody plant dominance, woody plant height (and both evergreen and deciduous heights), soil moisture, and soil organic carbon stock variables to reduce the heteroscedasticity in the component models. Before the log-transformation, woody plant dominance is inverted due to negative skewness, and after the log-transformation, it is inverted back. The transformation decisions are based on the visual interpretation of the histograms of the linear model residuals for each response variable. Due to the high Spearman correlation (− 0.72) between the snow depth and TPI (Supplementary material 2), we exclude TPI from models that use snow depth as a predictor (that is, all vegetation and soil models).

Secondly, we fit all hypothesised paths (hereafter, component models) using linear models in the SEM. The paths include both the direct pathways from the predictors (topography) and mediators (snow depth and vegetation) to the response variables (soil properties) (Figure 1). The SEM has seven component models: (1) snow depth predicted by elevation and TPI, (2–4) vegetation predicted by elevation, radiation, TWI, and snow depth, and (5–7) Soil properties predicted by all predictors and mediators (for a detailed model structure, see Supplementary material 3).

The model is based on the assumption that we can explain the overall physical foundation affecting the soil properties so well (Kemppinen and others 2018), that the remaining variation is mainly due to the local plant–soil interactions, that is, vegetation affecting the soil. We use high-resolution topography data and fine-scale in situ measurements of snow and vegetation to control all responses in our model. We account for all possible residual correlations among the vegetation variables as well as among the soil variables, as we assume that correlations may exist between the error terms.

Lastly, we build separate SEMs for evergreen and deciduous woody plants. The evergreen and deciduous dominances are highly correlated (− 0.70) (Supplementary material 2); therefore, we do not use them in the same model, and consequently, we build two separate SEMs. These models have the exact same model structure as described above, except for the woody plant dominance and woody plant height variables, which are replaced with either the evergreen or deciduous variables (for a detailed model structure, see Supplementary material 4).

Results

The 171 study plots captured large environmental gradients in topography (Figure 2), woody plant dominance and height (Figure 3), and soil properties (Figure 4).

Figure 4
figure 4

Tundra soil properties. Soil moisture (A) and temperature (B) data from study plots (n = 171) shown as box plots in relation to the weather conditions. Daily precipitation shown as bars (A). Daily mean air temperature from two observation stations presented as lines, with their ranges shown as polygons (B). The dashed line shows 0 °C (B). Data on the soil organic carbon stocks (C) presented as a violin plot overlaid with a box plot. In the box plots (A–C), first, the notches and the hinges represent the 25th, 50th, and 75th percentiles. Second, the whiskers represent the 95% confidence interval, with outliers shown as points. Third, the thickness of the boxes indicates the duration (in days) of each soil moisture and temperature measurement campaign (A, B). In the violin plot (C), the thickness of the violin polygon corresponds to the local density of the observations. Daily air temperature and precipitation data from the Saana observation station (1; 1002 m asl) and Kyläkeskus observation station (2; 480 m asl) (FMI 2019a, b).

Full size image

The R2 values for the seven component models in the SEM ranged from 0.09 to 0.54 (Figure 5; for detailed results of the model, see Supplementary material 3). The highest R2 values of the component models were found for snow depth (0.54), woody plant dominance (0.42), and soil moisture (0.42). According to tests of directed separation, we did not exclude any significant (P < 0.05) paths from the model, except for the link from TPI to woody plant height.

Figure 5
figure 5

Woody plants mediating the effects of topography and snow on tundra soils. Results of a structural equation model with seven component models. Here, we only present relationships with P < 0.05. The arrows show the standardised regression slopes associated with the links. Double headed arrows represent correlations between error terms. For results on evergreen and deciduous models as well as further details and relationships with P ≥ 0.05, see Supplementary material 3 and 4.

Full size image

In the SEM, the strongest links were found in the snow and vegetation component models: from topography to both snow depth and vegetation (Figure 5; Supplementary material 3). TPI correlated negatively with snow depth (standardised coefficient = − 0.70; P < 0.001), elevation with woody plant height (− 0.47; < 0.001), and TWI with woody plant dominance (− 0.41; < 0.001), respectively. In addition, snow depth negatively correlated with plant coverage (− 0.22; P < 0.001) and woody plant dominance (− 0.35; < 0.001).

We found significant links from plant coverage and woody plant dominance to all three soil properties (Figs. 5, 6; Supplementary material 3). We found that the component model on soil moisture had three links: TWI (standardised coefficient 0.35; P < 0.001), plant coverage (0.18; 0.031), and woody plant dominance (− 0.39; < 0.001). In the component model on soil temperature, we found five links: elevation (− 0.39; < 0.001), radiation (0.37; < 0.001), TWI (− 0.17; 0.048), plant coverage (− 0.02; 0.024), and woody plant dominance (− 0.22; 0.029). We found that the component model on soil organic carbon stocks had two links: plant coverage (0.19; 0.046) and woody plant dominance (− 0.34; 0.002). However, we did not find links from snow depth or woody plant height in any of the component models regarding the soil properties.

Figure 6
figure 6

Relationships between tundra vegetation and soil properties based a structural equation model. The variables are log-transformed, except for soil temperature. The points represent the partial residuals plotted between the vegetation and soil properties. The line represents the linear relationship between the variables. The statistically significant relationships (P < 0.05) are presented in colour, whereas the insignificant ones in grey.

Full size image

In the SEM with only evergreen dominance and height, we found one link from evergreen woody plants to the soil properties: evergreen height to soil organic carbon stocks (standardised coefficient 0.26; 0.002). In the SEM with only deciduous dominance and height, we found one link from deciduous woody plants to the soil properties: deciduous dominance to soil moisture (− 0.23; 0.001). For detailed results of SEMs with either only evergreen or deciduous woody plants, see Supplementary material 4.

Discussion

We investigated with a hierarchical modelling approach based on SEM, if the coverage and height of woody plants affect the growing-season topsoil moisture and temperature despite the potentially strong influence of topography and snow. We used in situ measurements from 171 plots combined with high-resolution LiDAR data on topography. We hypothesised that in addition to the influence of topography, snow accumulation, and the overall vegetation cover on growing season soil moisture and temperature (< 10 cm depth) as well as on soil organic carbon stocks (< 80 cm), the effect of these drivers on soils would also be manifested (that is, mediated) via the dominance and height of woody plants (Figure 1). We assumed that the effects of dwarf shrubs would be subtle, due to the strong links from topography to snow and vegetation. However, we found mediating effects of woody plants, as the dominance of woody plants decreases soil moisture, soil temperature, and soil organic carbon stocks (Figs. 5, 6).

We also investigated the effect of woody plant height on the soil properties; however, the model did not find links between these variables. We must acknowledge the possibility that the height of particularly the deciduous woody plants may also be inhibited by grazing (Vowles and others 2017; Maliniemi and others 2018). However in our plots, the height of the two dominant deciduous woody plants Betula nana (mean height 11.9 cm, range 1–45 cm) and Vaccinium myrtillus (7.1 cm, 2–25 cm) (Figure 3) fall within their commonly described heights across the tundra biome (14.6 cm, 1–79 cm; 8.6 cm, 1–26 cm) (Bjorkman and others 2018).

We tested the hypothesis separately for the evergreen and deciduous species and found links from evergreen height to soil organic carbon stocks and from deciduous dominance to soil moisture (Supplementary material 4). In the Fennoscandian mountain tundra, woody plants are chiefly evergreen dwarf shrubs (< 50 cm tall) (Vowles and Björk 2019), which grow not only as dense mats, but also rather sporadically intertwined with deciduous species. We speculate that separating the effects results in weaker relationships between woody plants and soils, which in turn, breaks links found when modelling the groups together. However, the results of the separate models add valuable insights into the mechanisms explaining how different functional groups of dwarf shrubs influence soil properties.

Unexpectedly, we found no direct effects of wintertime snow depth on the growing season soil microclimate or the soil organic carbon stocks. Previous studies found that wintertime snow conditions affect multiple soil properties, particularly growing season soil moisture (Ayres and others 2010), even when controlling for the effects of topographical variation (Williams and others 2009). The reason for the lack of this relationship is likely related to the woody plants reflecting this information already due to the strong plant-snow interactions in tundra environments (Niittynen and others 2018).

Effects of woody plants on soil moisture

Our investigation captured a large gradient in soil moisture throughout the growing season (2.2–83.4 VWC%) (Figure 4), which was the soil property most accurately modelled in our study (Figure 5). The model found links to soil moisture from topography, both directly and indirectly through the mediating effect of the dominance of woody plants. Also, wintertime snow depth was linked to soil moisture, but only through woody plants. We found that the overall dominance of woody plants as well as the dominance of deciduous woody plants directly related to less soil moisture.

Our findings are in support of previous studies, which found that the expansion of woody plants may lead to an increased evapotranspiration and intercepting precipitation (Bonfils and others 2012; Zwieback and others 2019). Previous studies reached conclusions similar to our fine-scale investigation: that is, soil moisture is generally lower in woody plant habitats compared to other tundra habitats (Ge and others 2017; Lafleur and Humphreys 2018).

Evergreen woody plants, such as Empetrum nigrum, can form dense matts leading to less evaporation, since they cast shade on the soil beneath them throughout the growing season. Extensive woody plant expansion can increase water vapour in the atmosphere, because woody plants transpire water more efficiently than barren tundra (Swann and others 2010) and possibly also more than other vegetation types. An increased transpiration, in turn, may amplify the greenhouse effect of woody plants alongside changes in the albedo (Swann and others 2010).

In this hierarchical setting, we found a relationship from vascular plant coverage to soil moisture; however, the link was weak and positive, unlike the influence of woody plant dominance on soil moisture, which we found to be strong and negative (Figure 6). This demonstrates that we were able to (1) separate the effects of woody plants from the overall effects of vascular plants and (2) distinguish between the effects of woody plants on soil moisture from the effects of soil moisture on vascular plants.

Effects of Woody Plants on Soil Temperature

In the beginning of the growing season, soil temperature was highly variable across plots (0.1–15.0 °C) and homogenised towards the end of the growing-season (3.1–7.9 °C) (Figure 4). Soil temperature was strongly linked to elevation and solar radiation (Figure 5). In addition, the model suggests a mediating effect of the overall vascular plant coverage and the dominance of woody plants, which both decrease soil temperatures (Figure 6). This is in support of previous studies, which have found that in treeless ecosystems such as the tundra, vegetation in general cools the soil surface by shading (Blok and others 2010; Lantz and others 2013; Myers-Smith and Hik 2013).

Woody plants may form dense mats, which effectively insulate and shade the soil during the growing-season, and in turn, decrease soil temperatures (Blok and others 2010). However, the shading effect of deciduous woody plants is dependent on phenology (Bonfils and others 2012), as they may reach their maximal shading potential only after their buds have turned into leaves. In addition, weather conditions may potentially affect the influence of dwarf shrubs on soil microclimate. We speculate that both phenology and weather may partly result in the shrunken soil temperature gradient towards the end of the growing season.

Woody plants can have strong effects on radiation interception and the mixing of air (Aalto and others 2018); thus, abundant vegetation may sustain a microclimate decoupled from the macroclimate (Asmus and others 2018). In the tundra, tall, erect species, such as Betula nana or Juniperus communis (Figure 3) can potentially cause wind turbulence, and, in turn, enable the release of latent heat from the soil to the open atmosphere decreasing the soil temperatures (Bonfils and others 2012). Whereas small, prostrate species such as cushion plants grow closer to ground and are more aerodynamic compared to taller species, which may lead to less turbulence at the ground surface. Thus, the effects of woody plants on soil temperatures are explained by several mechanisms associated with shading and mixing of air.

Effects of Woody Plants on Soil Organic Carbon Stocks

Soil organic carbon stocks were linked to topography, but only through the mediating effects of vegetation (Figure 5). We found no other links to the soil organic carbon stocks, which indicates that the model may have lacked other important factors explaining the stocks. However, we found a positive link from the overall plant coverage and the height of evergreen woody plants to soil organic carbon stocks and a negative link from the dominance of woody plants (Figure 5).

Our results are in support of previous literature, which have found that particularly deciduous shrubs produce litter that induces soil priming (see for example, Parker and others 2015; Sørensen and others 2018). Evergreen shrubs may increase soil organic carbon stocks (Vowles and Björk 2019), for example, by producing slowly decomposing leaf litter and recalcitrant organic complexes (Cornelissen 1996). Our results are also in support of findings, in which deciduous and evergreen shrubs growing in a warmer microclimate were net CO2 sources during the growing season due to their high respiration compared to herbaceous vegetation that was a CO2 sink (Cahoon and others 2012), suggesting that the carbon stocks of woody plants might be decreasing.

Our results can be further explained by the differences between woody and herbaceous plants. For example, plant root properties and rhizosphere processes can influence soil organic carbon stocks (Sistla and others 2013; Parker and others 2015). Root longevity (that is, age) is lower in graminoids than in shrubs (Iversen and others 2015), which can result in slower carbon cycling and lower carbon inputs in shrub understory soils. In addition, woody plant roots are shallower than graminoid roots; thus, woody plant expansion may result in smaller carbon inputs from roots to the deeper soil layers, and in turn, decrease the stocks (Ylänne and others 2018). However, understanding these various mechanisms would require a broader range of measurements beyond this study, particularly on soil organisms and roots.

Conclusion

We investigated the effects of woody plants on multiple soil properties in the tundra in a multivariate setting. Our results suggest that the dominance of woody plants decreases soil moisture, soil temperature, and soil organic carbon stocks. The structural equation modelling approach enables taking into account the hierarchy of this multivariate setting, considering both topography and wintertime snow conditions and their effects on vegetation and soil properties. In this study, we provide field evidence on the impacts of woody plants on tundra soils to (1) support findings in previous studies and (2) encourage future researchers to consider SEM for controlling the hierarchy within a study system. Considering the results presented here and in previous literature, as well as the magnitude of the expected expansion, woody plants and their expansion are likely to impact the entire tundra system through the water, energy, and carbon cycles of tundra. Thus, it is important to investigate woody plant expansion in relation to soil microclimate and carbon stocks.