Material and methods

Study area and plot layout

The study was performed across the eastern and western slopes of the Blæsedalen valley outside the town of Qeqertarsuaq, Disko Island, Western Greenland. The landscape features a heterogeneous topography, basaltic soils, and dwarf shrub-dominated tundra vegetation of variable density. We selected a total of 90 locations stratified randomly, aiming to cover gradients of elevation, vegetation productivity and water availability, across an area of 15 km2 (von Oppen et al. 2022; Figure 2a). At plot establishment, we buried a soil temperature logger (HOBO MX-2201, Onset Corp., Bourne, MA, USA) in 6 cm depth at 10 cm into true South direction from the plot centre. Additionally, we placed a TMS-4 microclimate logger (TOMST, Prague, Czech Republic; Wild et al. 2019) at 10 cm into true North direction from the plot centre, which recorded topsoil moisture content (Figure 2b). We established the plots between 26 June and 01 July 2019 and retrieved tea bags and soil temperature loggers between 19 and 27 August 2020, with the length of incubation period varying by a maximum of 12 days across plots.

Figure 2:

Study design and environmental variation. a) Study area on Disko Island, Western Greenland. White circles indicate plot locations and numbers show elevation strata used for stratified random plot placement (separated by black contours). b) Positioning of loggers (TMS-T4: black dot; HOBO: black diamond) and teabags (triangles) within sample plots. Arrangement of tea types was randomised between plots. c) Variation of topographic (brown), vegetation/soil (green) and microclimate (blue) variables across sample plots. Horizontal jittering of points indicates the density of the value distribution. Grey points in the background for the gridded topographic variables show values for 1000 randomly sampled points within the study area (representing the full environmental variation in the region).

Topography, vegetation and soil variables

We recorded position and elevation at the centre of each plot location using a Trimble R8s dGNSS system (Trimble Inc., Sunnyvale, California, U.S.). We also measured local slope direction and inclination across a 4 m-diameter circle using a handheld inclinometer, from which we calculated local solar radiation as mean Solar Radiation Index (SRI; Keating et al. 2007) across the study period. In addition, we used the ArcticDEM (digital elevation model) at 2 m resolution (version 3.0; Porter et al. 2018) to extract (i) Topographic Position Index (TPI; Weiss 2001) within a 30 m radius, and (ii) Topographic Wetness Index (TWI; Beven and Kirkby 1979), following recommendations by Kopecký et al. (2021), through the SAGA GIS (Conrad et al. 2015) TPI and TWI modules, respectively, for our plot locations.
During the 2019 growing season, we visually quantified the proportional cover of shrubs, graminoids, forbs, bryophytes and lichens within a circle of 50 cm radius around the plot centre (Figure 2b) and subsequently summed their cover values into a single vegetation cover variable. Values for vegetation cover could surpass 100 % in case of multiple vegetation layers. In addition, we sampled functional traits for a subset of the most dominant species in the 60 plots on the eastern slope. As dominant species we considered those with >5 % cover in one or more of the first 18 plots surveyed across low, mid and high elevation levels, plus additional locally abundant species with >5 % cover in any given plot. For each species within each plot, we selected up to three (if available) healthy individuals closest to the plot centre and within the plot radius. For these individuals, we measured both vegetative and maximum plant height as the vertical distance from the ground to the tallest vegetative plant part or tallest plant part irrespective of tissue type. Furthermore, we sampled up to three leaves, or a small branch for small-leaved dwarf shrubs, from each individual for leaf trait assessment (n = 522 samples). We then determined the leaves’ fresh weight and scanned them for leaf area calculation through ImageJ (Schindelin et al. 2012), before drying leaves at 70 °C for 48 hours. Using the leaves’ dry weight, we calculated Specific Leaf Area (SLA) as dry weight per area, and Leaf Dry Matter Content as leaf dry weight divided by fresh weight. Due to limited resources, we could only determine chemical traits for a subset of 300 samples, with one random sample per species, elevation band and approximate water availability category (see sampling design in von Oppen et al. 2022). Dry soil and leaf samples were milled. We then used 5 to 20 mg per sample, depending upon organic matter content, to measure total C and N as well as isotopic ratios of C (13C/12C) and N (15N/14N). These analyses were carried out using an elemental analyzer by Dumas combustion at 1020 °C (Flash 2000, Thermo Scientific, Bremen, Germany) coupled to a Thermo Delta V Advantage IRMS in continuous flow.
We imputed community-level trait moments for the remaining 30 plots using hierarchical nonparametric bootstrapping through thetraitstrap package (Maitner et al. 2021) incorporating the plots’ moisture and productivity classes as determined during stratified plot placement. We subsequently performed a principal component analysis (PCA) from which we extracted plot scores along the first and second main axis of trait variation, which were related to leaf economics (PC1) and plant size (PC2), respectively (Figure S1). These axes aligned well with previous findings for global (Díaz et al. 2016) and tundra plant functional trait space (Thomas et al. 2020).
We took soil samples 10 cm north (“central”) as well as at 50 and 150 cm distance into all intercardinal directions from the plot centre using a metal corer with 6 cm diameter and 3 cm depth. Samples were only taken if there was enough soil to fill at least 50 % of corer volume. Samples were oven-dried at 45 °C to constant mass and sieved to obtain gravel (> 2 mm diameter) and soil (< 2 mm) fractions. If less than 40 g of dry soil were available for a central sample, we selected one of the outer samples with sufficient material (n 50 cm = 30, n 150 cm = 16). If that was not the case for any of the nine samples, we pooled soil across samples, recording the contribution from each location (n = 22). Subsamples of 10 g were used to determine contents of several macro- (carbon, nitrogen, phosphorus) and micronutrients (Ca, K, Mg, Na), as well as soil pH and cation exchange capacity. Particle size fractions were determined from the remaining 30 g of soil (see Table S1 for a list of instruments used for soil analyses). We then performed a principal component analysis including all soil variables and, for each plot, extracted the scores for the first and second axis of variation. The two axes mainly corresponded to nutrient content, and soil pH and texture (higher silt, lower sand content), respectively (Figure S2).

Soil temperature and moisture data

We set soil temperature loggers to record soil temperature at 15-minute intervals throughout the tea bag incubation period. After visual assessment of temperature and moisture curves, we removed measurements from known periods of disturbance (field observations) during the 2019 growing season (± 1 day). These occurred for example if animals had removed loggers from the ground (n = 5). We also excluded records entirely if temperature loggers had been removed from the ground between 2019 and 2020 (n = 5). For each plot, we eventually calculated temperature sums (growing-degree days above 0 °C, GDD0) across the incubation period.
Soil moisture loggers recorded at 10-minute intervals over the 2019 growing season (27 June – 07 August). We removed values from known periods of logger disturbance (n = 2) and where extremely low moisture signals indicated insufficient contact of the probe with the soil matrix during the measurement period (n = 2). We calibrated soil moisture using soil type-specific calibration curves provided by the manufacturer, assigning soil types to plots based on minimum Euclidean distances for particle size fractions. As we obtained a few negative moisture values for some plots due to imperfect calibration, we scaled soil moisture values to the range of individual measurements. We then calculated mean soil moisture across the 2019 growing season for each plot.

Burial and processing of tea bags, calculation of decomposition metrics

To measure decomposition dynamics, we used commercially available Green Tea and Rooibos Tea (EAN 8714100770542 and 8711327514348; Lipton, Unilever Inc.), on the basis of the TBI protocol (Keuskamp et al. 2013). The TBI has been used previously in numerous locations across the Arctic tundra (Sarneel et al. 2020, Gallois et al. 2022, Björnsdóttir et al. 2022, Thomas et al. 2023). The method should therefore be well suited to investigate the effects of multiple landscape-scale environmental gradients on litter decomposition.
We labelled and dried tea bags at 70 °C for 48 hours before recording the initial dry weight of each individual bag. During plot establishment, we buried tea bags at 2 cm depth and 10 cm distance into intercardinal directions from loggers (i.e., NW and NE of microclimate loggers and SW and SE of temperature loggers; Figure 2b), in order NW - NE - SW - SE irrespective of tea type to randomise positioning, corresponding to two replicates per type per plot and a total of 360 tea bags. Exact position and/or depth could vary slightly in response to local soil conditions. We retrieved tea bags after approx. 14 months (416 to 428 days). Due to occasional damage or disturbance by animals, we collected 285 intact bags (81 %). The number of plots with at least one intact sample was 83 for Green Tea, 79 for Rooibos Tea, and 79 for both Green Tea and Rooibos Tea. We dried the bags at 70 °C for 72 hours and placed them in ziplock bags with silica gel for transport and storage. Starting in August 2021, we dried tea bags again at 60 °C for 24 hours to remove any moisture potentially drawn during the storage period. We then determined post-burial dry weight and subtracted the average weight of 10 empty bags from all weights to obtain initial and remaining weight as well as mass loss of tea for each sample. Finally, we calculated relative mass loss for Green Tea and Rooibos Tea for each sample, and mean values per plot.
The TBI (Keuskamp et al. 2013; see therein for all parameters and equations in this subsection) describes litter decomposition dynamics according to a two-phase model - fast decomposition during the initial phase and slow decomposition during the second phase. Green Tea as a labile substrate is characterised by a high hydrolysable fraction (\(H_{g}\) = 0.842). Green Tea therefore has a relatively short first phase of decomposition, which enables estimation of how quickly a labile substrate stabilises in a given environment. The stabilisation factorS quantifies this amount of undecomposed labile litter, with high values of S indicating less favourable conditions for decomposition. It is calculated as:
\begin{equation} S=1-a_{g}\cdot H_{g}^{-1}\nonumber \\ \end{equation}
with \(a_{g}\) as relative mass loss of Green Tea. Second, Rooibos Tea as a more recalcitrant material has a lower hydrolysable fraction (\(H_{r}\) = 0.552) and its first phase of decomposition therefore stretches over a longer time period. However, assuming that S is constant for both substrate types, the theoretically decomposable fraction of Rooibos Tea (\(a_{r}\)) can be calculated as:
\begin{equation} a_{r}=H_{r}\cdot(1-S)\nonumber \\ \end{equation}
This can then be used to calculate the decomposition rate constantk , which describes how fast this gradual decay process occurred for the decomposable litter fraction:
\begin{equation} k=ln(\frac{a_{r}}{M_{r}-(1-a_{r})})\cdot\frac{1}{t}\nonumber \\ \end{equation}
Here, \(M_{r}\) is the relative remaining mass of Rooibos Tea, and \(t\)the length of the incubation period in days. High values of ktherefore indicate rapid turnover of litter.

Tundra-wide decomposition data

To put decomposition rates from our study area into context, we obtained data on mass loss, k and S from the TundraTea initiative, which encompasses coordinated studies applying the TBI protocol across Arctic and alpine tundra regions (Thomas et al. 2023). We only included data from year-round incubations in untreated plots comparable to our study, and from regions with at least three unique plots to ensure variation in environmental conditions. This resulted in a total of 191 plots from 11 regions (1318 individual tea bags; Table S2).

Statistical analyses

From the 79 plots with complete decomposition response values, we removed: one extreme outlier with 2.4 times the next largest value fork , one plot which had no trait values due to sparse vegetation cover, and six plots with erroneous or missing microclimate measurements (resulting n = 71).
We used Structural Equation Modelling (SEM) to analyse the effects of variation in topography, vegetation and soil, and microclimate on decomposition dynamics. We chose SEM as the method is well suited to test our hypothesised system of direct and indirect processes influencing tundra soil decomposition, while also allowing us to account for any potential latent relationships not specified by us (Grace et al. 2015). According to our hypothesis (Figure 1), we fitted eleven sub-models (Table S3) describing the local impacts (i) of topography (elevation, solar radiation, topographic position and wetness) on vegetation (vascular plant cover and community-level leaf economic and size trait composition) as well as soil (nutrients and texture); (ii) of topography, vegetation and soil on microclimate (year-round GDD0 and growing-season soil water content); and (iii) of plant community traits, soil and microclimate on each of the decomposition response variables, i.e. mass loss of either tea type, as well as on decomposition rate k and stabilisation factorS . We then combined the sub-models into a single multiple linear regression model within the SEM, using the piecewiseSEM package (Lefcheck 2016) in R. We inspected residual distributions for all sub-models to ensure assumptions about homoscedasticity were met. We also calculated pairwise linear relationships between decomposition response and explanatory variables within sub-models for comparison (i.e., not incorporating any covariates, residual variation from other models, or effect size standardisation; Table S4). These individual relationships were consistent with SEM results in terms of magnitude and direction of relationships as well as significance level (Tables S3, S4), hence we describe and discuss only SEM results in the following, as they better represent our initial hypothesis.
All data management and analyses were performed in R v4.1.1 (R Core Team 2021).