Environmental drivers of increased ecosystem respiration in a warming tundra

Abstract

Arctic and alpine tundra ecosystems are large reservoirs of organic carbon^1,2. Climate warming may stimulate ecosystem respiration and release carbon into the atmosphere^3,4. The magnitude and persistency of this stimulation and the environmental mechanisms that drive its variation remain uncertain^5-7. This hampers the accuracy of global land carbon-climate feedback projections^7,8. Here we synthesize 136 datasets from 56 open-top chamber in situ warming experiments located at 28 arctic and alpine tundra sites which have been running for less than 1 year up to 25 years. We show that a mean rise of 1.4 °C [confidence interval (CI) 0.9-2.0 °C] in air and 0.4 °C [CI 0.2-0.7 °C] in soil temperature results in an increase in growing season ecosystem respiration by 30% [CI 22-38%] (n = 136). Our findings indicate that the stimulation of ecosystem respiration was due to increases in both plant-related and microbial respiration (n = 9) and continued for at least 25 years (n = 136). The magnitude of the warming effects on respiration was driven by variation in warming-induced changes in local soil conditions, that is, changes in total nitrogen concentration and pH and by context-dependent spatial variation in these conditions, in particular total nitrogen concentration and the carbon:nitrogen ratio. Tundra sites with stronger nitrogen limitations and sites in which warming had stimulated plant and microbial nutrient turnover seemed particularly sensitive in their respiration response to warming. The results highlight the importance of local soil conditions and warming-induced changes therein for future climatic impacts on respiration.

Fig. 1. Study area, duration and environmental context of warming experiments across the tundra biome.

a–c, ER was collected in 136 growing season measurement years across 56 independent OTC warming experiments, covering 28 distinct geographical sites across the tundra biome (a), capturing a wide range of experimental warming durations at the time of ER measurements (b) and environmental heterogeneity (c). The radius of the dots in a indicates the number of years in which ER was measured per experiment (range 1–13 years; 63% of experiments with range greater than 1 year), the colour of the dots indicates the number of observations per dataset, that is, number of ER measurements in the growing season across all years in which measurements were conducted (details in Extended Data Fig. 1, Supplementary Methods 1 and Supplementary Tables 1 and 2). The distribution of datasets is shown across the duration of experimental warming before ER measurements were taken, grouped into four classes (duration = ER measurement year − OTC treatment start year; note that ‘0’ refers to the first summer of experimental warming) (b) and categorical variables describing environmental contexts (c), from left to right: zone; soil moisture class; soil pH class; and vegetation type (B, barren; G, graminoid; P, prostrate shrub; S, erect shrub; W, wetland tundra). The maps were made with R.

Fig. 2. Effects of experimental OTC warming on ecosystem ER.

Experimental warming increased ER across the tundra biome but the magnitude of the response varied across time and space. Effect of OTC warming on ER Hedges’ SMD calculated as (mean ER of the warmed plots − mean ER of the control plots)/pooled standard deviation across the 136 growing season datasets (that is, unique experiment × ER measurement year combinations). On the top of the graph, a blue diamond shows the mean estimate (est. = 0.57 and 95% CI [0.44–0.70], error bars) of the ER response across the 136 datasets, as well as the Q value testing for heterogeneity and P value from the meta-analysis. Black dots represent ER Hedges’ SMDs of individual datasets and 95% CIs (black error bars) in alphabetical and chronological order. Individual datasets are represented by the experiment ID in black (left) and ER measurement year (right) in a colour scale ranging from dark blue, light blue, orange to red which represents increasingly longer warming duration at the time of ER measurements. Experiments with more than 1 year of ER data are grouped. See Supplementary Tables 1, 2 and 4 for details on the datasets and SMD and CI values. The black dashed vertical line (SMD = 0) represents no change in ER with warming whereas the areas to the right and left of it represent increased (SMD > 0) versus decreased (SMD < 0) ER with warming. Dashed vertical lines (x = 0.2, 0.5, 0.8 and −0.2, −0.5, −0.8) reflect small, medium and large positive and negative Hedges’ SMD, that is, increasingly greater ER increases and decreases with warming.

Fig. 3. Temporal patterns in ER responses with experimental warming duration.

ER responses to warming faltered around the end of the first decade of warming but did not wane in the long term. a, Mean ER response to warming (ER Hedges‘ SMD) across the four experimental duration age classes, showing the kernel density estimation of the underlying distributions in violin plots as well as the single-factor metaregression model estimates and 95% CIs (in black). b, Relation between experimental warming duration and ER response to warming across the four age classes (Extended Data Table 2) with 95% CI (light grey area). Hedges’ SMDs for individual datasets (grey bubbles) are shown in a and b, calculated as (mean of the warmed plots − mean of the control plots)/pooled standard deviation. Bubble size denotes the weight of the observation used in the metaregression model, quantified as the inverse of the square root of within-study variance, with greater bubbles indicating greater weights. Qm values, P values and sample sizes (n) are shown for the regression models, with Qm representing the Q value of importance of duration age class across all data in a and of experimental warming duration per age class in b. Significant regression lines in b are shown in blue. The black dashed horizontal line (SMD = 0) represents no change in ER with warming whereas the areas above and below represent increased (SMD > 0) versus decreased (SMD < 0) ER with warming. Dashed horizontal lines (y = 0.2, 0.5, 0.8 and −0.2, −0.5, −0.8) reflect small, medium and large positive and negative Hedges’ SMD, that is, increasingly greater ER increases and decreases with warming. For detailed model output, see Extended Data Table 2.

Fig. 4. Indirect warming effects on the ER response to warming.

a,b, Warming-induced changes in TN concentration (a) and soil pH of the mineral layer (b) affected the ER response to experimental warming. Warming-induced changes in TN and pH on the x axis are quantified as Hedges’ SMDs of the soil conditions between the warmed and control plots, for example, (mean TN of the warmed plots − mean TN of the control plots)/pooled standard deviation. ER Hedges’ SMDs for individual datasets are shown on the y axis with grey bubbles, calculated as (mean ER of the warmed plots − mean of the control plots)/pooled standard deviation. Bubble size denotes the weight of the observation used in the metaregression quantified as the inverse of the square root of within-study variance, with larger bubbles indicating greater weights. The significant regression lines with 95% CI are shown with blue lines and shaded areas, respectively. Bottom left in each panel shows the Qm (Q value of importance of the environmental drivers), P value of the metaregression model and the sample size (n, number of datasets). The black dashed horizontal line (y = 0) represents no change in ER with warming whereas the areas above and below represent increased (SMD > 0) versus decreased (SMD < 0) ER with warming. Dashed horizontal lines (y = 0.2, 0.5, 0.8 and −0.2, −0.5, −0.8) reflect small, medium and large positive and negative Hedges’ SMD, respectively, or increasingly greater ER increases and decreases with warming. The black dashed vertical line (x = 0) reflects no change in the soil condition with warming (with areas right and left of it representing increased versus decreased conditions with warming). For detailed model output, see Extended Data Table 3.

Fig. 5. Context-dependencies in the ER response to warming.

a,b, TN concentration (%) (a) and C:N ratio (b) of the mineral soil layer affected the ER response to experimental warming. The environmental drivers on the x axis reflect mean values of measured soil conditions at the control plots of each experiment with available data. ER Hedges’ SMDs for individual datasets are shown on the y axis with grey bubbles, calculated as (mean ER of the warmed plots − mean of the control plots)/pooled standard deviation. Bubble size denotes the weight of the observation used in the metaregression quantified as the inverse of the square root of within-study variance, with larger bubbles indicating greater weights. The significant regression lines with 95% CI are shown with blue lines and blue shaded areas, respectively. Bottom left in each panel shows the Qm (Q value of importance of the environmental drivers), P value of the metaregression model and the sample size (n, number of datasets). The black dashed horizontal line (y = 0) represents no change in ER with warming whereas the areas above and below represent increased (SMD > 0) versus decreased (SMD < 0) ER with warming. Dashed horizontal lines (y = 0.2, 0.5, 0.8 and −0.2, −0.5, −0.8) reflect small, medium and large positive and negative Hedges’ SMD, respectively, or increasingly greater ER increases and decreases with warming. For detailed model output, see Extended Data Table 4.

Fig. 6. Spatial patterns and uncertainty in the sensitivity of ER to warming across the tundra.

Maps showing spatial extrapolation of metaregression results for warming effects on ecosystem respiration (ER) and associated uncertainties across the arctic and circumarctic alpine tundra. a, Spatial distribution of predicted relative changes in ER resulting from 1.4 °C air warming. Values were obtained by combining the significant metaregression multifactor model with global gridded soil data for TN concentration and C:N ratio of the mineral layer (Methods). Plotted values are predicted mean values expressed as percentage change from present ER levels. For the whole region, ER increased from 3.4 to 4.3 PgC yr⁻¹ (increase of 0.86 PgC yr⁻¹ with s.d. of 1.36 PgC yr⁻¹) or by 25% (s.d. = 40%). b, Uncertainty of estimated ER changes expressed as coefficients of variation (standard deviations of model predictions divided by mean predicted values). These standard deviations are computed with Monte Carlo error propagation incorporating uncertainty associated with (1) metaregression model parameters and (2) soil database input values (Methods). c, Relative contribution of model parameter uncertainty versus input soil data uncertainty to Monte Carlo estimates of standard deviation. Error propagation was conducted twice, first combining both uncertainty sources as in b and subsequently by fixing metaregression parameters to their mean estimates and only allowing soil input data to vary. The plotted values show the ratio of the resulting standard deviations, that is, soil input uncertainty only/combined uncertainty. Colours closer to green indicate grid cells for which uncertainty in ER responses to warming are primarily driven by error in model metaregression model parameter estimates; colours closer to purple indicate grid cells for which ER prediction uncertainty is driven more by imprecision in soil database data. The maps were generated using Natural Earth (https://www.naturalearthdata.com/).

Extended Data Fig. 1. Maps showing the locations of the 28 sites (red dots) used for meta-analysis across four different continents: America (top left), Europe (top right), Asia (bottom left) and Australia (bottom right).

Details in Supplementary Table 1. The radius of the red dots reflects the number of years that ecosystem respiration (ER) data was measured at each site, resulting in a total of 136 datasets.

Extended Data Fig. 2. Uniform distribution of ER response to warming across the environmental context of all datasets.

ER response to experimental warming across the spatial environmental context (See Fig. 1 in Main text). Violin plots of the actual data across the categorical environmental drivers, that is, Climate Zone (a), Soil moisture class (b), Soil pH class (c) and Vegetation class (d), are displayed, showing the kernel density estimation of the underlying distributions. ER Hedges SMDs for individual datasets are displayed with grey bubbles, calculated as (mean ER of the warmed plots - mean ER of the control plots)/pooled standard deviation. Bubble size denotes the weight of the observation used in the metaregression models quantified as the inverse of the square root of within-study variance, with greater bubbles indicating greater weights. Within the violin plots, single-factor metaregression model estimates and 95% confidence intervals are displayed with black circles and error bars. Above the x-axis, the ‘Qm’ value represents the importance of the moderator or environmental driver with p-value (‘p-val’) for each model. Number of datasets (‘N’) for each environmental driver per category is shown above the violin plots. The black dashed horizontal line (SMD = 0) represents no change in ER with warming while the areas above and below represent increased (SMD > 0) vs. decreased (SMD < 0) ER with warming. Dashed horizontal lines (y = 0.2, 0.5, 0.8 and −0.2, −0.5, −0.8) reflect small, medium and large positive and negative Hedges SMD effect sizes or increasingly greater ER increases and decreases with warming. For detailed model output, see Extended Data Table 4.

Extended Data Fig. 3. Warming-induced changes in soil conditions drive ER response.

Trends in warming-induced changes in soil conditions and in the vegetation community driving the ER response to experimental warming (p < 0.1): Hedges SMD of bulk density (BD) (a) and C:N ratio (b) of the soil mineral layer and of total N (TN) concentration of the soil organic layer (c) and aboveground biomass (d) and graminoid cover (e) of the vegetation community. ER Hedges SMDs for individual datasets are displayed on the y-axis with grey bubbles, calculated as (mean ER of the warmed plots - mean of the control plots)/pooled standard deviation. Bubble size denotes the weight of the observation used in the metaregression quantified as the inverse of the square root of within-study variance, with larger bubbles indicating greater weights. Top left in each panel shows the sample size (‘N’, number of datasets) and bottom left shows the ‘Qm’ (Q-value of importance of the environmental drivers) and ‘p’-value of the metaregression models. The black dashed horizontal line (y = 0) represents no change in ER with warming while the areas above and below represent increased (SMD > 0) vs. decreased (SMD < 0) ER with warming. Dashed horizontal lines (y = 0.2, 0.5, 0.8 and −0.2, −0.5, −0.8) reflect small, medium and large positive and negative Hedges SMD, respectively or increasingly greater ER increases and decreases with warming. The black dashed vertical line (x = 0) reflects no change in the environmental condition with warming (with areas right and left of it representing increased vs. decreased conditions with warming). For detailed model output, see Extended Data Table 3.