Environmental drivers of increased ecosystem respiration in a warming tundra

Data collection

Study design

We quantified and analysed the effect of warming on ecosystem CO₂ respiration (ER) across tundra ecosystems by collecting original, previously published and unpublished, data from OTC warming experiments in the tundra biome. OTCs are small greenhouses which passively increase air temperatures during the snow-free season while allowing free entry of precipitation⁴⁴. They are commonly used to simulate climate warming at a plot-scale in low-stature alpine and arctic tundra systems (International Tundra Experiment network, ITEX⁶⁸). We contacted potential data contributors (1) through email using the networks ITEX, WARM, InterAct and Permafrost Carbon, (2) at the 20th ITEX Meeting 9–13 September 2019 in Parma, Italy, and (3) by retrieving contact details of authors of earlier published meta-analyses of ER measured in warming experiments. ER datasets were only included if (1) experiments were located in alpine or arctic tundra, (2) original ER measurements were available from both OTC and unmanipulated control (ambient temperature) plots and (3) ER was measured during the plant growing season (that is, June through August for all sites, except for the Australian site AUS_1: October through February; Supplementary Fig. 2). This resulted in 24,035 data points from 136 ER datasets (that is, a unique site × experiment × ER measurement year; Supplementary Tables 1 and 2) from 56 unique OTC experiments across 28 tundra sites.

One site could consist of one or several warming experiments (for example, covering local differences in vegetation communities or soil conditions) and contribute ER measurements from one or several occasions during one or several years (Supplementary Table 1 and Supplementary Figs. 1 and 2). Experiments differed in abiotic (for example, permafrost presence, ambient climate and soil type) and/or biotic (for example, vegetation and size of microbial communities) conditions, as well as in warming treatment specifications. First, at the time of ER measurements, experiments differed in duration of their warming treatment (0–25 years, mean 7). Second, despite the common experimental warming setup, experiments differed in OTC size (0.3–0.7 m height), duration of OTC deployment (for example, year-round placement or winter removal of OTCs) and methodology of flux measurements. We therefore collected methodological data on (1) OTC height, (2) whether or not OTCs were removed during winter, (3) CO₂ analyser producer (Los Gatos Research, LICOR, PP-system or Vaisala), (4) CO₂ analyser type (infrared or laser), (5) measurement system (automated or manual and closed or open), (6) the timing of ER measurements (daytime or night time) and (7) ER measurement plot size. Details on the procedure used to check for methodological biases introduced by the factors described above are presented in the Supplementary Discussion 2. Because we observed no significant effects of these methodological factors on the ER response, we assume there is no methodological bias in our results.

Ecosystem respiration

For each dataset, we requested plot-specific daily average ER data and standardized all individual data points to the same unit, that is, gCO₂ m⁻² d⁻¹. We removed 339 outliers (0.8% of the data) to account for measurement errors (data points outside the range of mean ER ± 3 × s.d. in a dataset following ref. ⁶⁹). Subsequently, we calculated growing season average ER per experimental treatment (control versus warmed) per dataset as grams of carbon in CO₂ m⁻² d⁻¹. For details on the sampling frequency, see Supplementary Fig. 2 and Supplementary Table 2. All experiments measured daytime ER using dark or opaque chambers, except for ALA_1, which used automated clear chambers. For this experiment, we extracted night-time measurements based on the photosynthetically active radiation values (PAR < 5 µmol m⁻² s⁻¹) and included these as ER measurements.

Effect size calculation

We calculated effect sizes⁴⁶ of the mean growing season ER response to warming, separately for each dataset. We used Hedges’ g SMD as primary effect size across all models and the log-transformed⁵ ROM as additional effect size to quantify percentage change in ER with warming. Hedges‘ SMD for each dataset was calculated as the (mean growing season average ER in warmed plots − mean growing season average ER in control plots)/pooled standard deviation. This is deemed an appropriate effect size to use here because it overcomes data measured on different scales through standardization, as well as corrects for variance heterogeneity that may be introduced with small sample sizes (n < 20) by including the pooled standard deviation^46,70. Mean Hedges’ SMD values for each dataset were accompanied by 95% CIs of the standard errors, which reflect the precision of the effect size for each dataset⁷⁰. The ROM was calculated as the log-transformed mean growing season average ER in warmed plots/mean growing season average ER in control plots, which we transformed to a mean percentage increase as 100×(exp(ROM) − 1), that is, a more intuitive measure of changes in ER. See Supplementary Table 2 for the number of day-average ER measurements and number of plots used for each of the 136 datasets and see Fig. 2 and Supplementary Table 4 for an overview of mean Hedges’ SMD and 95% CIs per dataset.

Temporal patterns

To assess the ER response to warming over time, we quantified the duration of experimental warming for each dataset as the difference between the year of ER measurement and the year in which the warming treatment started (Supplementary Fig. 1 and Supplementary Table 1). We further categorized the datasets into four ‘age classes’ on the basis of warming duration, with class 1 including datasets in which 0–5 years of warming had been ongoing at the time of ER measurements and classes 2, 3 and 4 in which 5–10, 10–15 and more than or equal to 15 years of warming had been ongoing, respectively (Fig. 1). Note that duration = 0 refers to when less than 1 year of warming took place, that is, one single summer of warming.

Environmental drivers ER

To quantify environmental drivers that might explain the variation in the growing season ER response to warming, we obtained several types of environmental data on a site-, experiment- or ER-measurement level and used these to analyse indirect warming effects and context-dependencies as drivers of variability in ER response to warming. See Supplementary Table 5 for sample sizes of all environmental drivers used to assess indirect warming effects and context-dependencies and Supplementary Methods 1 for a detailed explanation on how the different drivers were obtained and calculated.

Indirect warming effects

To investigate whether and how changes in environmental conditions induced by the warming treatment influenced the ER response to warming, we quantified changes in the below-mentioned climatic and soil conditions as well as in vegetation and microbial communities for each dataset, by calculating effect sizes for these drivers: effect size calculation of ER data (above), for example, (mean air temperature warmed plots − mean air temperature in control plots)/pooled standard deviation. In doing so, we then tested whether changes in the environmental drivers due to warming influenced the changes in ER due to warming (see later section on ‘Metaregression’). The included environmental drivers for this part of our analyses are related to (1) climatic conditions at ER measurement time, that is, air temperature (°C), soil temperature (°C) and soil moisture (%); (2) soil conditions at plot level, that is, SOM (%), total C concentration (TC) (%), total N concentration (TN) (%), C:N ratio, bulk density (BD) (g cm⁻³) and pH from the mineral, and/or organic, layer; as well as organic layer depth (OLdepth) (cm); (3) the vegetation community at plot level, that is, percentage cover per functional group (graminoids, forbs, deciduous shrubs, evergreen shrubs, mosses and lichens), aboveground biomass and mean community height, at plot level; and (4) the microbial community at plot level, that is, proxies for bacterial and fungal biomass and derived fungal:bacterial (FB) ratios (Extended Data Table 3). Although ER can be strongly correlated to GPP of the vegetation, our database does not include GPP data, as its inclusion and comparison across time and across several experiments is challenging because of the light-dependency of photosynthesis and associated high spatiotemporal variability of GPP measurements. However, both processes are (partly) driven by the activity of the vegetation, which in turn is driven by direct responses to their environment (for example, microclimate) as well as by, for example, its biomass and community composition. We therefore argue that to increase our mechanistic understanding of the underlying ecological drivers of ER response to warming, it is more meaningful to incorporate the role of the response of plant growth to warming for ecosystem respiration by including (changes in) plant biomass and community composition, as well as NPP (see section ‘Context-dependencies’) as potential drivers of ER, as we do here. Note that we refer to the soil total C and N concentration throughout the manuscript, because these drivers were measured on solid (dried) soil with CNS element analysers. Therefore, they refer to the sum of concentrations of both inorganic and organic compounds of C and N, even though TC most likely reflects primarily the SOC concentration because tundra soils contain limited-to-no inorganic C.

Context-dependencies

To evaluate whether and how the environmental context influenced the ER response to warming, we included drivers quantifying the climate, soil conditions and vegetation community (Extended Data Table 4 and Supplementary Methods 1). For each dataset, (1) climate drivers included zone (low arctic, high arctic and alpine) and permafrost probability (values between 0 and 1) at site level and mean air temperature, soil temperature and soil moisture in the control plots at ER measurement time. Further, (2) soil conditions drivers comprised mean values of the SOM, TC, TN, CN, BD, pH and OLdepth data from the control plots per dataset (see above). In addition, we included soil pH class (pH low less than 5, medium 5–7, high more than 7), soil moisture class (dry, mesic and wet) and soil carbon stock (t ha⁻¹) here at the experiment level. Finally, (3) vegetation community drivers included a vegetation class (five categories: barrens (B), graminoid tundra (G, grasses and sedges), prostrate shrub tundra (P), erect shrub tundra (S) and wetlands (W), that is, based on the circumpolar arctic vegetation map (CAVM)⁷¹) and NPP (kgC m⁻² yr⁻¹) at the experiment level.

Statistical analyses

All analyses were performed using R v.4.0.3 unless specified otherwise⁷².

Modelling approach

Meta-analysis

To evaluate the effects of experimental warming on ER, we performed multivariate meta-analysis with the rma.mv function from the metafor R package⁷³, using the Hedges’ SMD of the growing season average ER data as primary effect size (Fig. 2 and Extended Data Table 1). In doing so, we obtained a mean pooled effect size of warming, taking into account the nestedness and repeated measurements in our 136 datasets. Each dataset was weighted by the inverse of its estimated sampling variance to increase both the precision of the estimated effect sizes and the statistical power in the meta-analyses⁷⁴. We also performed a meta-analysis on the ROM as secondary effect size (see previous section on ‘Effect size calculation’), which provided a mean pooled ROM of warming (Extended Data Table 1). Further, to test whether environmental conditions (climate, soil, vegetation and microbial: see previous section on ‘Environmental drivers’) significantly changed themselves because of the warming treatment, we performed separate meta-analyses testing the effect of warming on each of these drivers, using Hedges’ SMD as primary effect size and raw mean differences (that is, the difference in absolute values of the driver between treatment and control means) as secondary effect size because the latter provides estimates of absolute changes in the drivers (Supplementary Table 3).

Metaregression

Temporal patterns. To evaluate whether and how the ER response varied with duration of experimental warming, we performed two types of single-factor metaregression models. First, we tested whether duration (number of years of warming at time of ER measurements) influenced the magnitude of the ER response (ER Hedges’ SMD). These models were performed across experiments and both (1) across the four age classes of duration, to analyse overall long-term patterns of warming effects and (2) in each age class (0–4, 5–9, 10–14 and more than 15 years). This enabled to explore potential nonlinearity of patterns in ER response rates with warming over time (Fig. 3b and Extended Data Table 2a). Second, we tested (1) whether the age class influenced the ER response, that is, whether the mean ER response differed among the four age classes in sign or direction, as well as (2) whether the mean ER response per age class was significantly different from zero (Fig. 3a and Extended Data Table 2b). This enabled to evaluate the persistency (continued existence) of the ER response over time.

Indirect warming effects. To assess whether the variation in the ER response to warming was driven by indirect effects of warming on environmental conditions, we ran single-factor metaregression models using the Hedges’ SMD of the different environmental drivers (reflecting warming-induced changes in climate, soil, vegetation and microbial community conditions) as individual predictors and ER Hedges’ SMD as response (see previous section on ‘Indirect warming effects’; Fig. 4 and Extended Data Table 3).

Context-dependencies. To assess whether variation in the ER response to warming could be explained by context-specific environmental drivers related to climate, soil or vegetation, we performed single-factor metaregression models testing the influence of each of these drivers on ER Hedges’ SMD (see previous section on ‘Context-dependencies’; Fig. 5 and Extended Data Table 4).

We chose single-factor models to analyse the indirect warming effects and context-dependencies because the environmental driver data were unevenly spread across datasets (Supplementary Table 5) with different drivers missing from different datasets. Hence, multifactor models would have led to significant reductions in sample size and thus power. Before all metaregressions, we checked for confounding and/or collinearity issues through boxplots and correlograms.

Random-effects structure

Because we did not assume all datasets to have the same true effect size and to account for methodological differences across the datasets, we included a random-effects structure in all models⁵. Models were fitted with the following crossed and nested random effect structure: (~1 | experiment/dataset) + (~year | experiment). The first component (the multilevel structure) corresponds to a three-level meta-analysis, allowing the individual estimates of each observation (dataset) to be nested in the experiment grouping‐level^75,76. The second component (the autoregressive component) accounts for repeated samplings in the same site over several time points (years), by using a continuous autoregressive variance structure of true effects over time^76,77,78. We modelled such a structure by setting struct = ‘CAR’ in the rma.mv function from the metafor R package⁷³. This significantly improved the main meta-analysis model, that is, the model with the CAR structure had lower Akaike information criterion values than the model without; also, the acf-plot of computed ‘normalized or Cholesky residuals’ showed that autocorrelation was removed by adding the CAR. This random effect structure was then used for all meta-analyses and metaregression models.

Model visualization and validation

We considered and interpreted results as significant in our study if P < 0.05, as is commonly done in meta-analyses^79,80,81, whereas trends (P < 0.1) are described but not interpreted as significant (Extended Data Fig. 3). We visualized meta-analysis results through forest plots using the forest function from the metafor package⁸. Significant metaregression results were visualized by showing the actual observations with scatterplots and violin plots for continuous and categorical drivers, respectively; as well as the model predictions, that is, regression lines or group means ± 95% CIs, respectively. Models were validated by evaluating funnel and profile plots and testing the residuals for the assumption of normality and variance homogeneity^73,82. Funnel plots of the main meta-analysis showed asymmetry, indicating possible ‘publication bias’ but because we did not work with only published datasets (that is, not a meta-analysis sensu stricto) we argue that this detected asymmetry most likely originated from true heterogeneity rather than publication bias⁸³. We did not observe any other irregularities. Further, we calculated R² of the significant metaregression models using McFadden’s pseudo-R² by comparing log-likelihood values of the models to the one of the null model fitted through maximum-likelihood estimation. We report the omnibus test statistic Qm to estimate the amount of residual heterogeneity explained by each driver in the metaregression models, as well as the Q or QE value test statistic for heterogeneity for meta-analysis or metaregression models, respectively.

Supplementary analyses

To test the robustness of our results and to gain mechanistic insights which could help with further interpretation of our results, we performed several supplementary analyses (Supplementary Discussion). First, to answer whether the increased ER was driven by increased plant (autotrophic) or microbial (heterotrophic) respiration, or by both, we collected all the available data on ER partitioning in tundra warming experiments and performed meta-analyses on ER, autotrophic respiration and heterotrophic respiration for this subset (Extended Data Table 1b and Supplementary Discussion 1). Respiration partitioning measurements are destructive and hence alter ecosystem dynamics, making long-term measurements difficult. Hence, there is limited availability of continuous, simultaneous measurements of total ER and ER partitioning at the same site. Second, we analysed whether the differences in ER measurement methodology or OTC experiment setup resulted in methodological bias in our results (Supplementary Fig. 3 and Supplementary Discussion 2). Third, we analysed potential interactions of OTC microclimate effects on the ER response (Supplementary Discussion 3). Fourth, we assessed whether mismatching of ER with environmental driver data from different measurement years influenced our results through performing a sensitivity analysis, that is, rerunning our metaregression models on the basis of several restrictive sample size scenarios (Supplementary Fig. 4 and Supplementary Discussion 4). Fifth, we assessed whether the unbalanced sampling design from including two experiments with a much longer time series of repeated ER measurements (ALA_1 and GRE_6; Supplementary Table 2) than other experiments had a strong influence on the meta-analysis outcome; Supplementary Discussion 5). Finally, we analysed the temporal patterns in experiments as well, to see if they differed from our results when analysing these patterns across experiments (Supplementary Fig. 5 and Supplementary Discussion 6). Overall, these supplementary analyses confirmed the robustness of our results.

Spatial upscaling

We estimated spatial patterns in the sensitivity of ER to experimental warming and the resulting change in ER, by upscaling our findings on context-dependencies across the arctic⁶ and circumarctic alpine⁸⁴ tundra region (Supplementary Fig. 6). We performed the upscaling on the basis of the two context-dependent drivers which were significant predictors of ER responses to warming based on our single-factor metaregression models, that is, TN concentration and C:N ratio of the mineral soil layer. For the upscaling, a two-factor metaregression model was used (see details below). The detailed procedure is outlined below and visualized in Supplementary Methods 2.

Soil input data

As input for the upscaling, we downloaded the mean, 5th percentile and 95th percentile soil data (of TN, SOC concentration (%) and BD) for the study area with 250 m resolution from the ISRIC soil data hub⁸⁵. Given that the study area was about 10 million km², the 250 m resolution was not computationally feasible for Monte Carlo uncertainty analysis (see section ‘Monte Carlo uncertainty analysis’) and thus, we aggregated the data to 1 km resolution (mean over each 1 km grid cell, ignoring the grid cells without values). A supercomputer was used to perform the analysis on the basis of the aggregated 1 km resolution. We first extracted the mean TN and SOC concentration for the mineral layer of each grid cell by defining the mineral soil layers within 0–60 cm depth using a TN threshold of 0.01 g g⁻¹, a SOC threshold of 0.1 g g⁻¹ and a BD threshold of 1 g cm⁻³. We then calculated depth-weighted averages of each soil driver (TN and SOC) over the mineral layer for each 1 km grid cell. To estimate standard deviation (s.d.) for each grid cell, we used the 5th and 95th percentiles of the SOC and TN layers calculated as s.d. = (q.95 − q.05)/(2 × qnorm(0.95)), where q.95 and q.05 are 95th and 5th percentiles and qnorm is the inverse cumulative distribution function of a standard normal distribution. Because there was no dataset available for C:N ratio, we used the means and standard deviations of TN and SOC data to derive grid-cell-specific C:N ratios as SOC/TN, as advised by ISRIC soil data hub personnel. This derivation and the related uncertainty are described below in the section Monte Carlo uncertainty analysis.

Upscaling model

The regression model for upscaling was built with the log-transformed ROM as effect size (response variable), which allowed to calculate projected respiration as percentage change in ER with 1.4 °C warming (that is, the average warming achieved by the OTCs in our database), starting from a baseline ER. We used a two-factor metaregression model based on all ER datasets for which both TN concentration and C:N ratios of the mineral soil layer were available, with the following estimated coefficients: ROM = 0.05–0.16 × TN + 0.01 × C:N (QM = 6.7 (P < 0.05), QE = 274 (P < 0.001), n = 39).

Monte Carlo uncertainty analysis

To incorporate the uncertainty in the soil data and their relationship with ER changes into our upscaling estimates, we performed Monte Carlo simulations. To propagate uncertainty in the input soil data, we generated 100 values for TN and SOC for each 1 km² grid cell using the means and standard deviations derived from the ISRIC soil data. For each cell we sampled from a truncated multivariate random normal distribution (rtmvnorm from R package tmvtnorm), truncating the distribution of both variables to minimum of 5% of the mean value to avoid negatives. Because truncation alters the characteristics of the sample, we adjusted the means so that the mean of the truncated distribution matched the target mean obtained from ISRIC. The correlation between TN and SOC was set at 0.8273. For each sample of TN and SOC, a C:N ratio was calculated. We took the potential bias of using C:N ratio based on TN and SOC as a proxy for C:N ratio into account by fitting a linear regression between the reported C:N values from the database of this study and the calculated C:N values based on the reported TN and TC (Supplementary Fig. 8). We fitted 100 regressions to the scatter plot and used each regression to estimate the C:N ratio of each iteration (n = 100) in each grid cell.

To propagate more uncertainty derived from the estimated relationship between ROM and TN and C:N, we then used the 100 pairs of TN and C:N of each 1 km² grid cell as input for the upscaling model using predict.rma function from the metafor R package⁷³. This resulted in 100 predictions for ER percentage change (with associated prediction standard errors) for each grid cell. These 100 parameter sets were combined to produce an overall predicted value by taking the average, with the associated standard deviation computed by treating the prediction distribution as a mixture of the 100 prediction intervals. The standard deviation of the resulting distribution is given by:

where n is the number of simulated soil datasets and MEAN and SE are the vectors of the corresponding predicted means and standard errors.

This resulted in a predicted relative change in induced by 1.4 °C warming averaged over the 100 predictions pairs and the corresponding combined standard deviation, for each 1 × 1 km² grid cell (Fig. 6 and Supplementary Fig. 7). To estimate spatial patterns in the absolute change in ER, we then multiplied this relative change in respiration with spatially explicit baseline ER (obtained by summing gridded data of soil respiration⁸⁶ (data compiled from 1960 to 2011) and plant respiration⁸⁷ (data from 2015)).

We then estimated the relative contribution of model parameter uncertainty versus input soil data uncertainty to Monte Carlo estimates of the ER standard deviation (Fig. 6). Model predictions were conducted twice on the 100 input sets, first propagating both uncertainty sources (Fig. 6b) and subsequently by fixing metaregression parameters to their mean estimates and only allowing soil input data to vary. The same input data simulations were used in both runs to keep the data-related uncertainty constant. Figure 6c shows the ratio of the resulting standard deviations, that is, soil input uncertainty only/combined uncertainty.

Finally, we explored how well the distribution of our dataset compares to the distribution of gridded soil data used for upscaling. Thereto, we created a scatter plot for TN and C:N ratio in which we overlaid both the observed (from our dataset) and gridded values (gridded soil dataset used for upscaling) (Supplementary Fig. 9). Although the observed values from our dataset are slightly biased towards low TN values, overall, the observed values cover the main concentration of TN-C:N ratio distribution (about 64% when measured with Convex Hull polygon).