A systemic overreaction to years versus decades of warming in a subarctic grassland ecosystem.

Temperature governs most biotic processes, yet we know little about how warming affects whole ecosystems. Here we examined the responses of 128 components of a subarctic grassland to either 5-8 or >50 years of soil warming. Warming of >50 years drove the ecosystem to a new steady state possessing a distinct biotic composition and reduced species richness, biomass and soil organic matter. However, the warmed state was preceded by an overreaction to warming, which was related to organism physiology and was evident after 5-8 years. Ignoring this overreaction yielded errors of >100% for 83 variables when predicting their responses to a realistic warming scenario of 1 °C over 50 years, although some, including soil carbon content, remained stable after 5-8 years. This study challenges long-term ecosystem predictions made from short-term observations, and provides a framework for characterization of ecosystem responses to sustained climate change.

Ecosystems are mosaics of plants, animals and microorganisms that, when undisturbed, interact in a dynamic equilibrium (steady state) with each other and the physical environment connecting them[1,2]. Climate warming has been repeatedly shown to affect many biota and their interactions[3-5], often with consequences for ecosystem functioning[6,7]. Yet most observations of warming effects on ecosystems come from subsets of taxa[8,9], interactions[10] or processes[11-13], and we know little about the overall trajectory of ecosystems experiencing climate change. Indeed, notable exceptions to this rule have largely focussed on carbon cycle processes[12-14], and attempts to step beyond single metrics of ecosystems have relied on composite proxies (for example, community biomass[2,15,16]) or a priori decisions about desired baseline conditions[17], both of which ignore positive or negative covariance[12,18-20] among interacting components and risk distorting the view of an ecosystem’s steady state. No coherent framework exists to describe the interplay between the biotic and abiotic components of a warming ecosystem, leaving us unable to accurately forecast the future of ecosystems following decades to centuries of climate change.

An ecosystem’s relationship with temperature can take one of four forms, which can be conceptualised by resistance-resilience theory[20]. First, an ecosystem may be entirely resistant to temperature change. Second, an ecosystem may react rapidly and permanently to warming (low resistance, low resilience). Third, an ecosystem may initially resist warming, but be driven to an altered state by sustained or intense warming (high resistance, low resilience). Fourth, an ecosystem may react strongly to the onset of warming, but recover under prolonged warming (low resistance, high resilience). Both the duration and magnitude of warming will influence the nature of this relationship, in that warming of high intensity may result in a faster transition between an ecosystem’s ambient and warmed states. However, to date no observations have allowed direct comparison of how the numerous different components of an ecosystem react to temperature change, if at all, or in what sequence. This is compounded by the short lifetime of ecological experiments, most of which are less than 10 years old[12,13,21] and typically fall short of known lags and demographic processes in ecological systems[18,19]. There is thus a pressing need to assess warming effects on multiple components of an ecosystem collectively, and to examine their persistence over timescales relevant to Earth’s systems.

We used 128 measured variables representing a wide range of biotic (plants and soil organisms) and abiotic properties, pools and processes of a subarctic grassland to make a comprehensive, decadal-scale assessment of warming effects on an ecosystem. We exploited the longest known in situ warming experiment[22], which captures at least 50 years of soil warming (hereafter >50 years, long-term) using geothermal gradients (0 to 18 ºC above ambient temperature throughout the soil profile), coupled with similar geothermal gradients capturing 5-8 years of warming in the same landscape (hereafter 5-8 years, short-term). The large, stable, high-resolution temperature gradients and long warming duration offered by geothermal systems make them uniquely placed to give detailed mechanistic insight into the responses of ecosystems to sustained warming. In 2008, new geothermal gradients emerged in the same grassland as the long-term warmed gradients with similar ambient control plots, allowing us to compare the responses of the same ecosystem to 5-8 years versus >50 years of warming. Our approach was threefold. We first characterised how the ecosystem had reacted to >50 years of warming. We then determined whether the ecosystem showed the same response after 5-8 years of warming. Finally, we used these responses and associations between individual variables to construct a framework describing how warming affected the whole ecosystem.

Ecosystem response to >50 years of warming

We represented the ecosystem as the first axis (principal component; PC1) of an empirical orthogonal function (EOF) containing all observations and combinations of warming intensity and duration (see Methods). PC1 explained 33.7% of total variance, which was more than double the variance collectively explained by PC2 and PC3 (8.8% and 6.2%, respectively; see Supplementary Information). The variables with the highest loadings on PC1 described pools and processes throughout the plant-soil system, including the soil carbon stock, large water-stable soil aggregates, soil bacterial and fungal biomass, soil microbial community composition and plant stoichiometry, phenology and species richness (Supplementary Table S1). PC1 was thus a good representation of the ecosystem, but could not have been embodied by a single variable or several variables from the same subsystem (for example, aboveground biomass[2,15]). We found that PC1 was affected strongly by warming (LR = 68.87, df = 1,7, N = 59, P < 0.0001), but the nature of its response depended on warming duration (LR = 9.89, df = 1,9, N = 59, P = 0.0071). After >50 years, warming had a linear effect on the ecosystem (Fig. 1a). This held true for warming up to 18 ºC, encompassing and even widely surpassing the most severe climate scenarios for the next 300 years[23]. As such, the ecosystem did not resist sustained warming, but instead changed linearly with warming intensity, and no tipping points[24] in ecosystem structure or function were apparent with warming up to 18 ºC.

Figure 1

Whole-ecosystem responses to soil warming.

Responses of grasslands (N = 59) exposed to (a) long-term (>50 years, yellow) or (b) short-term (5-8 years, red) soil warming. Data are PC1 scores (33.7% explained variance) from a single empirical orthogonal function (EOF) containing 128 variables (see Methods). Statistics and fit lines reflect significance of warming (W), duration (D) and their interaction (W × D), as determined by GLS models (see Supplementary Table 2 for test outputs). (c) The reaction (Δ response) of the ecosystem to short-term warming, calculated as the difference between responses to short-term and long-term warming. Fit line is a loess smoothing function. In all panels, grey ribbons represent 95% confidence intervals of a null model testing for artefacts arising through data handling (see Methods).

Ecosystem response to 5-8 years of warming

We used the response to >50 years of warming as a benchmark model for contrasting short-term and long-term warming effects on the ecosystem. Ambient temperature plots were similar between short-term and long-term warmed grasslands (see Methods), and most (92 of 128) variables shared the same relationship with temperature after short-term and long-term warming (see below; Supplementary Table S1). As such, we considered the 5-8 year warmed ecosystem to represent an intermediate state that will converge on the long-term response after >50 years. Despite this, we found that short-term warming had a different (LR = 9.89, df = 1,9, N = 59, P = 0.0071), non-linear (Fig. 1b), effect on the ecosystem, whereby 5-8 years of warming up to 14 ºC had stronger effects on PC1 from the same EOF than >50 years of warming (Fig. 1c). This stronger short-term response was not a symptom of greater variability in the ecosystem’s initial relationship with warming, since coefficients of variation were consistent between short-term and long-term warmed transects when calculated either for PC scores or for variables individually (Supplementary Figs 4 & 5). These results demonstrate that the ecosystem did not resist warming on either timescale, and also that it overreacted to warming in the short-term. This overreaction was evident after 5-8 years of warming, but was lost entirely after >50 years. Thus, while the ecosystem remained different from control plots following >50 years of warming, it recovered partially from its initial reaction over 8 to 50 years. While temporary warming effects on components of ecosystems are not uncommon[1,25], we show here that overreactions to warming are systemic because they were detected in the most important axis of variation from a large set of ecosystem state and process parameters. Moreover, warming effects were sustained well beyond the lifetime of such an overreaction, not lost over periods of months or years (for example, ref. 25).

Grouped variables and their responses to warming

We grouped variables by their relationships with temperature to explore their individual roles in the overreaction from the ecosystem. Variables ranged from being unresponsive to warming (32% of variables, 16 plant-related, 8 microbe-related, 15 soil properties, 2 ecosystem fluxes; Extended Data Fig. 1) to showing one of three response types (temporally consistent, overreacting, under-reacting; Fig. 2), which we summarised using the first PCs of EOFs performed separately on each group. Considering positive (Fig. 2a-c) and negative (Fig. 2d-f) relationships together, 15% of variables (11 plant-related, 4 soil biota-related, 4 soil properties) responded more strongly to warming of 5-8 years than >50 years (Fig. 2b,e). Although these variables were only a subset of those measured, they were responsible for eliciting the same warming response from the EOF of the full ecosystem. These variables included ephemeral increases in microbial activity[1], plant phenology[26] and plant carbon to nitrogen ratios, temporary shifts in some aspects of soil fungal community composition[27] and attenuating losses of root, soil fungal and soil bacterial biomass (Fig. 3b,e). Thus, while the ecosystem as a unit overreacted to warming, its response was governed by a subset of components therein. Overreacting variables almost exclusively (15 of 19 variables) described biotic pools and processes, yet none are measures classically used to represent ecosystem stability or functioning[2,15,16] and most are absent from even the most comprehensive assessments of warming effects on ecosystems[12-14,28]. It is clear from this finding that an ecosystem’s biota play a fundamental role in dictating the temporal variability of its response to warming, but this role is often overlooked by a priori decisions about which components of ecosystems to track over time.

Extended Data Fig. 1

Unresponsive group under soil warming

The response of grouped variables exposed to long-term (>50 years, yellow; LT) or short-term (5-8 years, red; ST) warming (N = 59). Data are PC1 scores from an EOF performed on the group displaying no significant responses to sustained warming. Statistics reflect significance of warming (W), duration (D) and their interaction (W × D), as determined by a GLS model.

Figure 2

Response shapes under soil warming.

(a-c) Positive and (d-f) negative responses of grouped variables exposed to long-term (>50 years, yellow/dashed; LT) or short-term (5-8 years, red/solid; ST) warming (N = 59 in all cases). Data are PC1 scores from EOFs performed separately on groups displaying (a,d) stable (ST = LT), (b,e) overreacting (ST > LT) and (c,f) under-reacting (ST < LT) responses to warming (see Fig. 3 for individual responses). Statistics and fit lines reflect significance of warming (W), duration (D) and their interaction (W × D), as determined by GLS models (see Supplementary Table 2 for test outputs). Yellow/dashed and red/solid lines illustrate LT and ST responses, respectively, and black lines illustrating the response where no significant W × D interaction occurred. Inlays show reactions (Δ responses) to ST warming, calculated as for Fig. 1c.

Figure 3

Variable groups and their responses to warming.

Positive (left) and negative (right) responses of ecosystem properties, pools and processes to short-term (ST; 5-8 years) and long-term (LT; >50 years) warming (N = 20 in all cases). Variables were manually grouped by relationships with temperature (see Supplementary Table S1): (a,d) permanent/stable (ST=LT; orange), (b,e) overreactions (ST>LT; red), (c,f) underreactions/buffered (ST

We found that 39% of variables (17 plant-related, 9 soil biota-related, 24 soil properties, 1 ecosystem flux) displayed the same response type after 5-8 years and >50 years of warming (Fig. 2a,d). This group included rapid and temporally consistent shifts in soil abiotic properties, the composition of plant and soil microbial, and particularly bacterial, communities[27] and declines in the soil carbon stock and other organic matter pools[29] (Fig. 3a,d). These variables did not reflect the ecosystem’s overreaction to warming, but instead were stable after 5-8 years. Among this group were variables that described aspects of community composition and carbon cycling, raising confidence in existing assessments of warming-induced changes to ecosystem community structure[8] and soil carbon[12,13]. A final smaller group (13% of variables, 6 plant-related, 3 microbe-related, 8 soil properties) resisted 5-8 years, but not >50 years, of warming (Fig. 2c,f), and represented apparently buffered changes to some aspects of plant metabolism[30], stoichiometry and growth, alongside lagged declines in the richness of plant and soil fungal communities[27] (Fig. 3c,f). Short-term observations clearly underplay the influence of sustained warming on such variables. For example, we underestimated plant species losses by 6 to 11 species over 50 years of warming if we only used short-term data (see below). This prediction does not consider the arrival of novel plant species in the community, which may mask losses of extant species, but is more than the total species loss expected in similar ecosystems over the same timeframe[3].

A framework for the ecosystem’s response to warming

Our findings collectively suggest that the ecosystem had reached a new steady state after >50 years of warming, and that this steady state was dependent on warming intensity. This is because temporally consistent changes occurred to most components of the ecosystem (Fig. 2a,d) despite ephemeral (Fig. 2b,e) and delayed (Fig. 2c,f) changes to others, and the short-term response of the ecosystem matched its long-term response at high warming intensities (Fig. 1c). Both imply that the long-term response of the ecosystem is a temporally stable state that will not be surpassed by further warming, barring future evolution[31] or the arrival of new species in the community[10], and moreover that warming above 14 ºC will accelerate the convergence of the ecosystem to its warmed state.

With this, and the PC scores from the full ecosystem and groups of variables therein, we propose a sequence to the ecosystem’s response to warming. First, warming accelerates soil biotic activity[1] (for example, decomposition of litter and soil organic matter) and some aspects of plant physiology (for example, minimum NDVI), while also lengthening the growing season[26] (Fig. 3a,b,e). Accelerated biotic activity is facilitated by an abundance of soil carbon and nutrients, including previously stable pools of soil organic matter[29]. Most pools decline rapidly within 5-8 years (Fig. 3d), which changes the soil structure. Nevertheless, energy and nutrient exchange among biota, as well as declines in plant and soil microbial and nematode abundance (Fig. 3d,e)[1,32], create a transitory phase where elevated biotic activity persists after soil organic matter is depleted (here, still occurring after 5-8 years of warming). Such “ecological inertia” is temporary because it reflects a deficit between ecosystem supply and biotic demand, which selects against species with exploitative resource use strategies (for example, arbuscular mycorrhizal fungi; see [27]) and leads to community restructuring over 8 to 50 years (Fig. 3c,f). Our data suggest that long-term persistence under warming may be limited to species with the capacity to adjust their resource use strategies, such as some aspects of metabolism and elemental ratios (Fig. 3b,c,e,f). Ultimately, as the community changes, the ecological inertia is lost and biotic activity, while still accelerated, partially attenuates per unit of soil or area (Fig. 3b,e). The outcome for the ecosystem is the emergence of a new warmed state with a different soil and biotic composition that is again in balance with the biomass and activity of the biota therein (Fig. 1a)[1].

General implications

We have shown here that the outcome of warming for this ecosystem is characterised by both an initial reaction and its convergence to a less extreme long-term response. This is important because it suggests that ecosystem responses to warming may only become predictable after several decades, making inferences from short-term experiments challenging[12,13]. We confirmed this reasoning by testing whether the temporal dynamics we observed had a bearing over predictions of long-term ecosystem change under realistic rates of climate warming. Specifically, we calculated the potential magnitude of error generated when using short-term observations to predict the long-term responses of all measured variables to 1 ºC of warming, which corresponds to the magnitude of warming expected over 50 years under the IPCC’s most conservative climate change scenario (RCP 2.6). We found that short-term observations yielded predictions that were, on average, 124.6% larger than those arising from long-term observations (Fig. 4), translating to errors of greater than 50% for 113 out of 128 variables and errors of greater than 100% for 83 out of 128 variables. This exercise not only confirmed that large errors can be made when using short-term (here, 5-8 year) responses to make long-term predictions, but also revealed that even small responses to minor warming can have implications for an ecosystem when considered over timescales relevant to climate change. Crucially, without making such calculations we may have wrongly concluded that expected warming in this region will have a negligible effect on the ecosystem, given our observations that warming effects on PC scores became most evident with warming above 3 ºC (Figs 1 & 2). We posit that the apparent discrepancy between conclusions drawn from prediction errors (Fig. 4) and PC scores (Figs 1 & 2) arose due to heterogeneity in the ecosystem’s biota under ambient temperature conditions[33], leading to uncertainty regarding the ecosystem’s pre-warmed state. Indeed, we suggest that the large warming range exploited by this experiment helped to characterise the responses of variables to low intensity warming in the face of such heterogeneity, and to constrain resulting predictions over 50 years of expected climate change. Taken together, these findings provide evidence that warming effects on ecosystems are relevant at low warming intensities, irrespective of uncertainty around them or their associated statistical significance, and advocate consideration of timescales and temperature ranges that go beyond those captured by the majority of existing warming experiments.

Figure 4

Prediction errors from short-term observations.

The distribution of error (%) generated when making long-term predictions from short-term observations (N = 128). Error was calculated as the absolute discrepancy between long-term and short-term responses of all 128 variables to 1 ºC of warming, reflecting the change expected over 50 years under the most conservative IPCC climate scenario (RCP 2.6). The x-axis is on a log10 scale, with a value of 100 indicating a magnitude of error of 100%.

In summary, this study demonstrates a clear need to target indicators of both the temporal dynamics and future warmed state of an ecosystem to fully understand its response to temperature change. Variables related to soil microbial activity and plant phenology, which here overreacted to 5-8 years of warming[1,26], may be useful metrics for tracking an ecosystem’s trajectory following the onset of warming. At the same time, plant and soil community composition and the soil carbon stock, which here were stable after 5-8 years of warming, may be appropriate indicators of the likely state of an ecosystem experiencing sustained warming. Our results originate from a subarctic grassland exposed to two discrete timescales of warming, so we call for future work to interrogate these variables as potential proxies against the existing suite of warming experiments worldwide[12-14,28,33]. We also call for further studies to use such existing platforms to validate the sequence of the warming response we report here, perhaps in the context of a hierarchical response framework[34], with particular attention to how species richness in both plant and soil communities changes between 10 and 50 years. In conclusion, the framework presented here is the first timeline for simultaneously mapping many properties, pools and processes onto an ecosystem’s overall trajectory under temperature change. It also delivers a list of variables that separately describe the temporal dynamics and warmed state of an ecosystem experiencing long-term warming. We urge consideration of this framework in future assessments of climate warming impacts on ecosystem structure and functioning, including decadal- to centennial-scale feedbacks to Earth’s systems.

Methods

Site description

We made use of the geothermal warming sites of the ForHot experiment[22] near Hveragerdi in Iceland (64°00'01″ N, 21°11'09″ W, 83-168 m a.s.l.). The experiment is situated on unmanaged grasslands in two valleys dominated by Agrostis capillaris, Ranunculus acris and Equisetum pratense over a Brown Andosol of approximately pH 5.7. One valley has been warmed consistently for at least 50 years, but likely since records began in 1708[22] (>50 years; long-term), whereas the other has been warmed since an earthquake on 29 May 2008 (5-8 years; short-term). The valleys each contain five replicated soil warming gradients (50-100 m length) ranging from ambient temperature (mean annual soil temperature: 5 ºC) to + 20 ºC, all of which are associated with different geothermal sources (see Supplementary Fig. S9). Warming in all gradients is seasonally consistent and has been stable since measurements began in 201322. To avoid confounding effects of geothermal activity on soil hydrology, half of the gradients were established uphill from a heat source and the other half downhill from a heat source. No substances associated with geothermal activity have been found in any plot since the experiment began[1,22,27,29]. While short-term and long-term transects were situated in adjacent valleys that shared the same geology, climate and land use history, we caution that it is not possible to eliminate the potential for pre-existing differences between valleys to have influenced comparisons between them. Nevertheless, we found no evidence that such differences occurred, for four reasons. First, considering all 128 variables together, ambient temperature plots were as similar within the short-term and long-term warmed transects as between them (Euclidean Distances: LR = 1.18, df = 1,3, N = 45, P = 0.2765). Second, 122 out of 128 variables (95%) did not significantly differ between ambient temperature plots (Bonferroni-adjusted P > 0.05 in all cases, N = 10), with only soil small and large macro-aggregate contents, soil sulphur and aluminium concentrations and plant potassium and manganese concentrations differing in baseline conditions between short-term and long-term warmed transects. Third, PC1 scores from the full empirical orthogonal function (EOF) containing all variables and plots (see below) did not differ between the short-term and long-term warmed transects independently of warming intensity and prior to normalising baselines (LR = 0.51, df = 1,3, N = 59, P = 0.4742). Finally, PC1 scores from the EOF of the full ecosystem not only shared the same pre-warmed state, but also converged on the same state with warming above 14 ºC (Fig. 1). It is extremely unlikely that any pre-existing differences between valleys would be detectable under minor to moderate warming but be undetectable under ambient conditions or extreme warming. Given this, we considered ambient temperature plots to be equivalent across all transects irrespective of warming duration. While no experimental system is without limitations, our approach overcomes some major criticisms of warming experiments to date[12,13,21,35], specifically by considering two timescales of warming throughout the soil profile over a large warming range and in a regression-style design.

Data collection & pre-processing

We collected data representing the per plot relative abundances of 11424 soil bacterial/archaeal operational taxonomic units (OTUs), 1447 soil fungal OTUs, 16 soil microbial phospholipid fatty acid markers, 43 plant species and 52 metabolites from two plant species, as well as another 110 variables representing other properties, pools and processes of the plant and soil system (Supplementary Table S1). All measurements were taken between 2013 and 2016 and expressed according to standard protocols (see Supplementary Table S3). Response variables with more than 50% missing values were removed (19 variables). The four-year sampling period was small compared to the difference between 5-8 years and >50 years of warming, which was a minimum of 42 years. Nevertheless, climatic variation unrelated to the warming transects resulted in interannual variability in some multi-year measurements (see Supplementary Information). We accounted for this by measuring plant phenology, biomass and ecosystem CO2 fluxes, which are variables known to be seasonally variable, on multiple dates over the four-year period and expressing each as the plot-level mean of all dates. We also only considered variables collected for all plots within the same year(s). Ecosystem CO2 flux data were further corrected for unrelated covariance in abiotic variables by expressing them as the residuals of models including photosynthetically active radiation, soil moisture and excess soil temperature variation as explanatory variables (see Supplementary Information). We collapsed multivariate datasets, namely microbial community composition, plant community composition and Anthoxanthum odoratum and Ranunculus acris metabolism, to three axes of an ordination and a measure of richness each (see Supplementary Information). We standardised the final 128 variables by centring around the mean and dividing by two standard deviations, and expressed every variable as the within-grassland difference between plot values and the mean value of ambient temperature plots. This approach yielded three ecosystem states: (i) a non-warmed ecosystem; (ii) the ecosystem following 5-8 years of warming; and (iii) the ecosystem following >50 years of warming. Thus, we could characterize the temporal dynamics of warming effects on the ecosystem in a fully replicated design using plots possessing a numerically identical pre-warmed state. Finally, mean summer temperature (MST: May to September, 2013 to 2016) at 10 cm depth was derived from hourly records (HOBO TidbiT V2 Water Temperature Data Loggers; Onset Computer Corporation, USA) in each plot. We expressed warming (ºC) as the within-transect differences between a plot’s MST and the mean MST for ambient plots, and removed one plot with MST warming >20 ºC.

Representing the ecosystem

We expressed the full ecosystem as the first axis (PC) of an EOF containing a total of 128 variables. An EOF is functionally comparable to a Principal Component Analysis but is not constrained by the same assumptions and accepts missing values[36]. In doing so, we were able to consider a large number of state and process parameters simultaneously, with no a priori decisions about their weighting (with the exception of multivariate data, see Data collection & pre-processing, above) and irrespective of possible covariance among them[17]. The EOF yielded a similar ordination to a Principal Coordinates Analysis (PCoA; Extended Data Fig. 2), which is a classical ordination approach that also accepts missing values. PC1 of the EOF explained 33.7% of total variance, with PC2 and PC3 explaining 8.8% and 6.2%, respectively. Given the large decline in explained variance between the first and next axes, we considered PC1 scores to be a good representation of the ecosystem, but note that informative warming effects were also found on PC2 and PC3 (see Supplementary Information).

Extended Data Fig. 2

Comparison between an EOF and a PCoA

The effects of warming intensity (contours) and duration (colours, marginal boxplots) on a grassland exposed to long-term (>50 years) or short-term (5-8 years, red) soil warming (N = 59), as represented by the first two axes of an (a) empirical orthogonal function (EOF) and (b) Principal Coordinates Analysis (PCoA).

Grouping variables by their warming response

We explored the different types of response exhibited by components of the ecosystem by grouping variables based on their relationships with warming and summarising these groups using separate EOFs. Grouping was performed using a three-step process. First, we used the P-values of warming effects and warming × duration interactions (Supplementary Table S1) with an α cut-off of 0.05 to categorise variables as temporally dynamic (W × D: P < 0.05), temporally consistent (W × D: P > 0.05, W: P < 0.05) or unresponsive (W × D: P > 0.05, W: P > 0.05). Second, we used coefficients from the same models to attribute a positive or negative direction to warming effects. Finally, we visually inspected temporally dynamic variables to determine whether warming effects were larger or smaller in the short-term versus the long-term, which we described as overreactions or under-reactions, respectively. It is important to note that P-values were used to define a cut-off in the first step of this process, but were not used as evidence of statistically significant warming effects on individual variables. Rather, statistical significance was tested at the group level on PC1 scores from EOFs performed separately on each group (see Statistical analysis, below), and was interpreted with reference to the per ºC changes of individual parameters reported in Fig. 3 (see Plotting relationships, below). While grouping based on any criterion is subjective, we used P-values because they are an established metric for examining the probability of biologically meaningful relationships, are functionally equivalent to using likelihood ratios, effect sizes or coefficients and, as opposed to these alternatives, already possess accepted thresholds. Here, we considered an α threshold of 0.05 because histograms of P-value distributions for warming effects and warming × duration interactions showed that α = 0.05 marked a threshold below which the frequency of P-values increased (Supplementary Fig. 6b,d). Despite this, we performed a sensitivity analysis to determine how changing the α cut-off between 0.05 and 0.01 in 0.01 steps would alter the composition of groups (Supplementary Table S4). For variables categorised as temporally dynamic at α = 0.05, a change in the threshold to α = 0.01 resulted in 9 out of 36 moving to a temporally stable group and 4 out of 36 becoming non-responsive. For variables categorised as temporally stable at α = 0.05, a change in the threshold to α = 0.01 resulted in 15 out of 76 becoming non-responsive. Despite some reshuffling of variables between groups, PC1 scores from the EOFs performed on separate groups remained numerically similar at all α cut-offs tested (Pearson Product Moment correlation: r > 0.8 and P < 0.0001 in all cases; Supplementary Fig. S7), meaning that the impact on the grouping process was negligible and general patterns of over- and under-reactions held true irrespective of the cut-off chosen. We thus proceeded with the groupings arising from the most inclusive α threshold of 0.05, but caution that the membership of a variable to a particular group is not definitively proven and there is a risk of misclassification for a small number of variables (Supplementary Table S4). This approach yielded a total of seven groups representing positive and negative temporally consistent, overreacting and under-reacting responses to warming, in addition to the unresponsive group.

Statistical analysis

We used generalized least squares (GLS) models to determine the effects of warming intensity (ºC above ambient MST), warming duration (transects warmed for 5-8 or >50 years) and their interaction on the PC1 scores from the EOF representing the full ecosystem, as well as on the PC1 scores from the EOFs representing different groups of variables therein (Supplementary Table S2). We also used GLS models with the same structure to generate P-values for effects of warming intensity, warming duration and their interaction on variables individually (Supplementary Table S1), but note that statistics performed on individual variables were used only to assign variables to groups (see Grouping variables by warming response, above) and no corrections were made for multiple testing. While the frequency of significant P-values observed for warming and warming × duration effects was higher than expected based on chance alone (Supplementary Fig. S6a,c), we advise against interpreting effects on individual parameters without considering further P-value corrections. GLS models were used so that, where necessary, we could account for unequal variance in explanatory variables. We scrutinized GLS model fits using residuals versus fitted values plots, histograms of residuals and boxplots of residuals against individual explanatory variables. In all cases, models included warming as a second-order polynomial, which was simplified to a linear term if it was non-significant (P > 0.05). Test statistics were obtained using sequential single-term deletions followed by likelihood ratio tests between models including and excluding explanatory terms.

Plotting relationships

We plotted PC1 scores from all EOFs against warming intensity grouped by warming duration. We also derived a new variable to illustrate whether PC1 scores changed more or less under short-term versus long-term warming. Given that exact temperatures differed between all plots, this was achieved using fitted GLS models to predict values for a long-term response using the temperatures from the short-term warmed plots. We then calculated the difference between the short-term response and the expected long-term response, creating a variable representing the reaction to short-term warming. Positive reaction values indicated a larger response to 5-8 years than >50 years of warming, negative reaction values indicated a smaller response to 5-8 years than >50 years of warming, and reaction values of zero indicated no change in the response to warming between 5-8 and >50 years. Finally, we expressed the changes of individual variables under warming using standardised (Fig. 3, graphs) and original-unit (Fig. 3, values in columns) per ºC changes over the full warming intensity range.

Validating relationships with null models

We performed a series of empirically-derived simulation analyses based on null models to rule out the possibility that observed effects on PC scores could have been caused by bias introduced through the data handling process. This was necessary for three reasons. First, original data was centred around the within-grassland means of ambient temperature plots. Using mean values alone ignored possible variance in the ambient temperature treatment, which has the potential to yield error in centred values that could amplify or dampen observed differences between warmed and ambient plots. Second, error in the centring process could have been incorporated into the EOF and, in an unlikely worst-case scenario, become the most important axis of variation (PC1) in ordinated data. While linear models formally include a null hypothesis that no relationship between X and Y exists, in light of such potential error it is not necessarily intuitive how PC scores would behave given no relationships with warming intensity or duration. Finally, it is not obvious how error and uncertainty surrounding real relationships with warming and PC scores would together proliferate through calculations of the ecosystem’s reaction to short-term warming, which we derived from observed and predicted PC scores. We thus used a simulation analysis to create null models based on 4000 randomised permutations that accounted for these potential sources of error. This was achieved by first calculating the within-grassland upper and lower 95% confidence intervals for ambient treatment mean values and creating four datasets that were centred around each of the four possible combinations of these intervals. Second, for each dataset separately we performed an EOF and created 1000 randomised permutations of the resulting PC scores. This yielded a total of 4000 sets of randomised PC scores that incorporated error arising from the centring process and for which we expected no relationships with warming intensity or duration. Third, we performed GLS models including warming intensity, warming duration and their interaction to predict values for all sets of randomised PC scores in 1 ºC steps over a 0 to 20 ºC warming range and in each grassland separately. We interpreted these predicted values as a set of 4000 null datasets for the effects of warming intensity and duration on PC scores. We then used randomised PC scores and corresponding predicted values from the same models to calculate the differences between short-term and long-term warmed grassland PC scores across a warming range of 0 to 20 ºC (see Plotting relationships, above). We interpreted these differences as a set of 4000 null datasets for the reaction of the ecosystem to short-term warming. Finally, we expressed both sets of null models as the 95% confidence intervals of all permutations therein, calculated separately for every 1 ºC warming step and, in the case of warming effects on PC scores, also separately for grasslands. In doing so, we used the largest possible confidence interval range for each warming step and grassland combination. We visualised 95% confidence intervals of null models as ribbons underlying corresponding figure panels for PC scores of the full ecosystem (Fig. 1) and an additional figure for PC scores of grouped variables (Supplementary Fig. S8). Overall, this process yielded null models that were in all cases unrelated to warming intensity or duration, eliminating the possibility that warming effects on PC scores were an artefact of the data handling process.

Prediction errors

We estimated the potential magnitude of error generated when using short-term observations to predict the responses of all variables to 1 ºC of warming. We selected 1 ºC because it reflected the magnitude of warming expected under the most conservative IPCC climate change scenario over a 50-year period (RCP 2.6), which is the minimum possible warming duration captured by the long-term warmed transects. We used predicted values from GLS models (see Statistical Analysis, above) to estimate the change in each variable between 0 and 1 ºC of warming separately for the short-term and long-term warmed transects. We then used these values to derive a prediction error (%) associated with the short-term responses, which we defined as the absolute percent difference between the changes in short-term versus long-term transects:

where ST and LT represent the predicted values for a variable at 1 ºC warming in the short-term and long-term warmed transects, respectively, and ST and LT are the predicted values for the same variable at 0 ºC in the short-term and long-term warmed transects, respectively. In doing so, we approximated the magnitude of error generated when making long-term predictions using data from short-term observations, which we plotted using a histogram and boxplot on a log10 scale (Fig. 4).

Supplementary Material

Supplementary Information is linked to the online version of the paper at www.nature.com/nature.