This is the accepted version of the article: Walker, Tom W. N.; Janssens, Ivan A.; Weedon, James T.; [et al.]. «A systemic overreaction to years versus decades of warming in a subarctic grassland ecosys- tem». Nature ecology & evolution, Vol. 4, Issue 1 (January 2020), p. 101–108. DOI 10.1038/s41559-019-1055-3 This version is avaible at https://ddd.uab.cat/record/232201 under the terms of the license A systemic overreaction to years versus decades of warming in a subarctic grassland ecosystem Tom W. N. Walker 1,2 *, Ivan A. Janssens 3, James T. Weedon 4, Bjarni D. Sigurdsson 5, Andreas Richter 6,7, Josep Peñuelas 8,9 , Niki I. W. Leblans3,5, Michael Bahn10, Mireia Bartrons 8,11, Cindy De Jonge 3, Lucia Fuchslueger 3,6 , Albert Gargallo-Garriga8,9,12, Gunnhildur E. Gunnarsdóttir5,13, Sara Marañón-Jiménez 3,8,9, Edda S. Oddsdóttir 14 , Ivika Ostonen 15 , Christopher Poeplau 16, Judith Prommer 6, Dajana Radujković 3, Jordi Sardans8,9, Páll Sigurðsson5, Jennifer L. Soong3,17, Sara Vicca 3, Håkan Wallander 18, Krassimira Ilieva-Makulec 19 and Erik Verbruggen 3 Abstract 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 Main 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 them1,2. Climate warming has repeatedly been shown to affect many biota and their interactions 3,4,5, often with consequences for ecosystem functioning6,7. Nevertheless, most observations of warming effects on ecosystems are derived from subsets of taxa 8,9, interactions10 or processes11,12,13 and we know little about the overall trajectory of ecosystems experiencing climate change. Indeed, notable exceptions to this rule have largely focused on carbon cycle processes12,13,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 covariance12,18,19,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 conceptualized 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 than10 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–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 provide 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 versus >50 years of warming. Our approach was threefold. We first characterized 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. Results and discussion 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). Variables with the highest loadings on PC1 described pools and processes throughout the plant–soil system, including 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 1). 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, above-ground biomass2,15). We found that PC1 was affected strongly by warming (likelihood ratio = 68.87, d.f. = 1,7, N = 59, P < 0.0001), but the nature of its response depended on warming duration (likelihood ratio = 9.89, d.f. = 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 of up to 18 °C, encompassing and even widely surpassing the most severe climate scenarios for the next 300 years (ref. 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. Ecosystem response to 5–8 years of warming We used the response to >50 years of warming as a benchmark model for contrasting short- and long-term warming effects on the ecosystem. Ambient temperature plots were similar between short- and long-term warmed grasslands (see Methods), and most (92 of 128) variables shared the same relationship with temperature after short- and long-term warming (see below; Supplementary Table 1). 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 (likelihood ratio = 9.89, d.f. = 1,9, N = 59, P = 0.0071), nonlinear (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- and long- term warmed transects when calculated for either PC scores or variables individually (Supplementary Figs. 4 and 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–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, underreacting; Fig. 2), which we summarized 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 phenology26 and plant carbon/nitrogen ratios, temporary shifts in certain 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 functioning2,15,16 and most are absent from even the most comprehensive assessments of warming effects on ecosystems 12,13,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. We found that 39% of variables (17 plant-related, nine soil biota-related, 24 soil properties, one 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, in particular, bacterial—communities 27 and declines in the soil carbon stock and other organic matter pools29 (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 carbon12,13. A final, smaller group (13% of variables, six plant-related, three microbe-related, eight soil properties) resisted 5–8 years—but not >50 years—of warming (Fig. 2c,f), and represented apparently buffered changes to certain aspects of plant metabolism 30, stoichiometry and growth, alongside lagged declines in the richness of plant and soil fungal communities27 (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 used only 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 exceeds the total species loss expected in similar ecosystems over the same timeframe3. 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 community10, and moreover that warming >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 certain aspects of plant physiology (for example, minimum normalized difference vegetation index, 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 ref. 27) and leads to community restructuring over 8–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 certain 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 characterized 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 become predictable only after several decades, making inferences from short-term experiments challenging12,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 Intergovernmental Panel on Climate Change’s (IPCC’s) most conservative climate change scenario (representative concentration pathway 2.6 (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 >50% for 113 out of 128 variables and errors of >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 >3 °C (Figs. 1 and 2). We posit that the apparent discrepancy between conclusions drawn from prediction errors (Fig. 4) and PC scores (Figs. 1 and 2) arose due to heterogeneity in the ecosystem’s biota under ambient temperature conditions33, leading to uncertainty regarding the ecosystem’s pre-warmed state. Indeed, we suggest that the large warming range exploited by this experiment helped to characterize 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. 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 probable 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,13,14,28,33. We also call for further studies to use existing platforms to validate the sequence of the warming response we report here, perhaps in the context of a hierarchical response framework34, 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 facilitates the simultaneous mapping of 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 the 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 above sea level). 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 probably since records began in 1708 (ref. 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. 9). Warming in all gradients is seasonally consistent and has been stable since measurements began in 2013 (ref. 22). 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- and long-term transects were situated in adjacent valleys sharing 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 among 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- and long-term warmed transects as between them (Euclidean distances: likelihood ratio = 1.18, d.f. = 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 sulfur and aluminium concentrations and plant potassium and manganese concentrations differing in baseline conditions between short- and long-term warmed transects. Third, PC1 scores from the full EOF containing all variables and plots (see below) did not differ between the short- and long-term warmed transects independently of warming intensity and before normalization of baselines (likelihood ratio = 0.51, d.f. = 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 >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 and preprocessing We collected data representing the per-plot relative abundances of 11,424 soil bacterial/archaeal operational taxonomic units (OTUs), 1,447 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 1). All measurements were taken between 2013 and 2016 and are expressed according to standard protocols (see Supplementary Table 3). Response variables with >50% missing values were removed (19 variables). The 4-year sampling period was small compared to the difference between 5–8 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 CO 2 fluxes, which are variables known to be seasonally variable, on multiple dates over the 4-year period and expressed each as the plot-level mean of all dates. We also considered only variables collected for all plots within the same year(s). Ecosystem CO 2 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 R. acris metabolism, to three axes of an ordination and a measure of richness each (see Supplementary Information). We standardized 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: (1) a non-warmed ecosystem, (2) the ecosystem following 5–8 years of warming and (3) 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–September, 2013–2016) at 10- cm depth was derived from hourly records (HOBO TidbiT V2 Water Temperature Data Loggers, Onset Computer Corporation) 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 PC 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 and preprocessing, above) and irrespective of possible covariance among them17. The EOF yielded a similar ordination to a principal coordinates analysis (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 following 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). 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 summarizing these groups using separate EOFs. Grouping was performed using a three-step process. First, we used the P values of warming effects and warming (W) × duration (D) interactions (Supplementary Table 1) with an α cut-off of 0.05 to categorize 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 over- or underreactions, 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 1 º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 W and W × D 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 steps of 0.01 would alter the composition of groups (Supplementary Table 4). For variables categorized as temporally dynamic at α = 0.05, a change in the threshold to α = 0.01 resulted in nine out of 36 moving to a temporally stable group and four out of 36 becoming non-responsive. For variables categorized 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, P < 0.0001 in all cases; Supplementary Fig. 7), meaning that the impact on the grouping process was negligible and that general patterns of over- and underreactions 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 that there is a risk of misclassification for a small number of variables (Supplementary Table 4). This approach yielded a total of seven groups representing positive and negative temporally consistent, overreacting and underreacting 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 2). We also used GLS models with the same structure to generate P values for effects of W, D and W × D on variables individually (Supplementary Table 1), but note that statistics performed on individual variables were used only to assign variables to groups (see Grouped variables and their response to warming, above) and no corrections were made for multiple testing. While the frequency of significant P values observed for W and W × D effects was higher than expected based on chance alone (Supplementary Fig. 6a,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- versus long-term warming. Given that exact temperatures differed among 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 than >50 years of warming, negative reaction values indicated a smaller response to 5–8 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 in individual variables under warming using standardized (Fig. 3, graphs) and original-unit (Fig. 3, values in columns) changes per 1 °C 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 were 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 4,000 randomized 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 we separately performed an EOF and created 1,000 randomized permutations of the resulting PC scores. This yielded a total of 4,000 sets of randomized 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 steps of 1 °C over a warming range of 0–20 °C and in each grassland separately. We interpreted these predicted values as a set of 4,000 null datasets for the effects of warming intensity and duration on PC scores. We then used randomized PC scores and corresponding predicted values from the same models to calculate the differences between short- and long-term warmed grassland PC scores across a warming range of 0–20 °C (see Plotting relationships, above). We interpreted these differences as a set of 4,000 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 visualized 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. 8). 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- 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 percentage difference between the changes in short- versus long-term transects: where STw and LTw represent the predicted values for a variable at 1 °C warming in the short- and long-term warmed transects, respectively, and ST a and LTa are the predicted values for the same variable at 0 °C warming in the short- 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 log 10 scale (Fig. 4). Reporting Summary Further information on research design is available in the Nature Research Reporting Summary linked to this article. Data availability Raw sequences (FASTQ format) are accessible through the NCBI Sequence Read Archive, under accession nos. SRP099121 and SRP075563 for bacteria (16S) and fungi (ITS1), respectively. Other data supporting the findings of this study are available in Figshare with the data https://doi.org/10.6084/m9.figshare.9958931 Acknowledgements T.W.N.W. is supported by the Swiss National Science Foundation (grant no. 31003A- 176044). J.L.S. is supported by the Office of Biological and Environmental Research, the US Department of Energy (contract no. DE-AC02-05CH11231). This work was further supported by a European Research Council Synergy Grant (no. ERC-2013- SyG 610028-IMBALANCE-P, awarded to I.A.J. and J.P.); a joint FWO–FWF grant (nos. FWO-G0F2217N and FWF-I-3237, awarded to I.A.J. and M.Bahn); three European Union Marie Skłodowska-Curie grants (nos. COFUND-291780 and Fellowship- 676108, awarded to S.M.J., and Fellowship-707270, awarded to C.d.J.); a Flanders Research Foundation Aspiration Grant (no. 11G1613N, awarded to N.I.W.L.); the Research Fund of the University of Antwerp (TOP-BOF and Methusalem grants, awarded to I.A.J.); grants from the Spanish Government (no. CGL2016-79835-P), the Catalan Government (no. SGR 2017-1005) and the Institut d’Estudis Catalans (no. PRO2008-SO2-PENUELAS) awarded to J.P.; the Icelandic Research Fund (no. 163272- 053 FORHOT-FOREST, awarded to B.D.S.); a JPI Climate Project (COUP-Austria, no. BMWFW-6.020/0008, awarded to A.R.); and the European Regional Development Fund (Estonia, Centre of Excellence ENVIRON and EcolChange). We also thank the Agricultural University of Iceland and Icelandic Forest Research for logistical support. References 1.Walker, T. W. N. et al. Microbial temperature sensitivity and biomass change explain soil carbon loss with warming. Nat. Clim. Change 8, 885–889 (2018). 2.Pennekamp, F. et al. Biodiversity increases and decreases ecosystem stability. Nature 563, 109–112 (2018). 3.Thomas, C. D. et al. Extinction risk from climate change. Nature 427, 145–148 (2004). 4.Walther, G.-R. et al. Ecological responses to recent climate change. Ecol. Lett. 416, 389–395 (2002). 5.Bardgett, R. D., Manning, P., Morriën, E. & De Vries, F. T. Hierarchical responses of plant–soil interactions to climate change: consequences for the global carbon cycle. J. Ecol. 101, 334–343 (2013). 6.Bragazza, L., Parisod, J., Buttler, A. & Bardgett, R. D. Biogeochemical plant–soil microbe feedback in response to climate warming in peatlands. Nat. Clim. Change 3, 273–277 (2012). 7.Giardina, C. P., Litton, C. M., Crow, S. E. & Asner, G. P. Warming-related increases in soil CO2 effux are explained by increased below-ground carbon flux. Nat. Clim. Change 4, 822–827 (2014). 8.Pearson, R. G. et al. Shifts in Arctic vegetation and associated feedbacks under climate change. Nat. Clim. Change 3, 673–677 (2013). 9.Blankinship, J. C., Niklaus, P. A. & Hungate, B. A. A meta-analysis of responses of soil biota to global change. Oecologia 165, 553–565 (2011). 10.Alexander, J. M., Diez, J. M. & Levine, J. M. Novel competitors shape species’ responses to climate change. Nature 525, 515–518 (2015). 11. Peñuelas, J., Rutishauser, T. & Filella, I. Phenology feedbacks on climate change. Science 324, 887–888 (2009). 12.Melillo, J. M. et al. Long-term pattern and magnitude of soil carbon feedback to the climate system in a warming world. Science 358, 101–105 (2017). 13.Crowther, T. W. et al. Quantifying global soil carbon losses in response to warming. Nature 104, 104–108 (2016). 14.Hicks Pries, C. E., Castanha, C., Porras, R. & Torn, M. S. The whole-soil carbon flux in response to warming. Science 355, 1420–1423 (2017). 15.Tilman, D., Reich, P. B. & Knops, J. M. H. Biodiversity and ecosystem stability in a decade-long grassland experiment. Nature 441, 629–632 (2006). 16.Hu, Z. et al. Shifts in the dynamics of productivity signal ecosystem state transitions at the biome-scale. Ecol. Lett. 21, 1457–1466 (2018). 17.Maestre, F. T. et al. Plant species richness and ecosystem multifunctionality in global drylands. Science 335, 214–218 (2012). 18.Bertrand, R. et al. Changes in plant community composition lag behind climate warming in lowland forests. Nature 479, 517–520 (2011). 19.Toth, L. T., Kuffner, I. B., Stathakopoulos, A. & Shinn, E. A. A 3,000-year lag between the geological and ecological shutdown of Florida’s coral reefs. Glob. Change Biol. 24, 5471–5483 (2018). 20.de Vries, F. T. et al. Land use alters the resistance and resilience of soil food webs to drought. Nat. Clim. Change 2, 276–280 (2012). 21.Wolkovich, E. M. et al. Warming experiments underpredict plant phenological responses to climate change. Ecol. Lett. 485, 21–24 (2012). 22.Sigurdsson, B. D. et al. Geothermal ecosystems as natural climate change experiments: the FORHOT research site in Iceland as a case study. Iceland Agric. Sci. 29, 53–71 (2016). 23.IPCC Climate Change 2013: The Physical Science Basis (eds Stocker, T. F. et al.) (Cambridge Univ. Press, Cambridge, 2013). 24.Scheffer, M., Carpenter, S., Foley, J. A., Folke, C. & Walker, B. Catastrophic shifts in ecosystems. Nature 413, 591–596 (2001). 25.Luo, Y. Q., Wan, S. Q., Hui, D. F. & Wallace, L. L. Acclimatization of soil respiration to warming in a tall grass prairie. Ecol. Lett. 413, 622–625 (2001). 26.Leblans, N. I. W. et al. Phenological responses of Icelandic subarctic grasslands to short-term and long-term natural soil warming. Glob. Change Biol. 23, 4932–4945 (2017). 27.Radujkovic, D. et al. Prolonged exposure does not increase soil microbial community response to warming along geothermal gradients. FEMS Microbiol. Ecol. 94, fix174 (2017). 28.Elmendorf, S. C. et al. Global assessment of experimental climate warming on tundra vegetation: heterogeneity over space and time. Ecol. Lett. 15, 164–175 (2012). 29.Poeplau, C., Kätterer, T., Leblans, N. I. W. & Sigurdsson, B. D. Sensitivity of soil carbon fractions and their specific stabilization mechanisms to extreme soil warming in a subarctic grassland. Glob. Change Biol. 23, 1316–1327 (2017). 30.Gargallo-Garriga, A. et al. Impact of soil warming on the plant metabolome of Icelandic grasslands. Metabolites 7, 44–22 (2017). 31.Walker, T. W. N. et al. Plastic and genetic responses of a common sedge to warming have contrasting effects on carbon cycle processes. Ecol. Lett. 22, 159– 169 (2018). 32.Marañón-Jiménez, S. et al. Geothermally warmed soils reveal persistent increases in the respiratory costs of soil microbes contributing to substantial C losses. Biogeochemistry 138, 245–260 (2018). 33.Langley, J. A. et al. Ambient changes exceed treatment effects on plant species abundance in global change experiments. Glob. Change Biol. 24, 5668–5679 (2018). 34. Smith, M. D., Knapp, A. K. & Collins, S. L. A framework for assessing ecosystem dynamics in response to chronic resource alterations induced by global change. Ecology 90, 3279–3289 (2009). 35.Carey, J. C. et al. Temperature response of soil respiration largely unaltered with experimental warming. Proc. Natl Acad. Sci. USA 113, 13797–13802 (2016). 36. Taylor, M. H., Losch, M., Wenzel, M. & Schröter, J. On the sensitivity of field reconstruction and prediction using empirical orthogonal functions derived from gappy data. J. Clim. 26, 9194–9205 (2013).