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).