Introduction

Northern terrestrial ecosystems face the fastest pace of warming and are considered among the most sensitive to climate change1,2,3,4. However, their responses are difficult to predict due to the highly variable sensitivity of their biophysical components and potential interactions among them. Some components can respond rapidly to warming (e.g., arctic lakes can shift to a new state within a few years5,6), while others are likely to react more slowly (e.g., vegetation shifting from tundra to forest-tundra over decades7,8,9). A lack of coherence in response time between components and across latitudes can restructure ecosystems and transform landscapes by modifying the combination of states among the coexisting components. To anticipate the impact of global warming on northern landscapes we must better assess and integrate the heterogeneity between components responding on varying timescales and their interactions.

The warming-induced state change of one ecosystem component can have abrupt and far-reaching impacts on other components. For instance, thawing permafrost can release massive amounts of carbon and nutrients, modify soil conditions for primary production, and affect geomorphological processes, hydrology, peatland dynamics, and human infrastructure10,11,12,13. Similarly, a shift in vegetation state from polar to shrub tundra can radically transform terrestrial ecosystems, offering new habitats and refuges for vertebrates and affecting microclimate and snowpack properties14,15. A change (or lack of change) in one component can trigger (or prevent) important changes in another component. Such cascading effects strongly influence the heterogeneity in response times and hence must be fully integrated in ecosystem sensitivity analyses (Fig. 1).

Quantifying the local and latitudinal heterogeneity in the response times of interacting ecosystem components poses a notable challenge for traditional scientific approaches. The complex cascade of events spreading at relatively large spatial scales prevents experimental and reductionist studies. Furthermore, approaches that rely on massive datasets to assess sensitivity can currently be applied at large spatial scales for very few components of the ecosystems1,16,17. Finally, cascading effects between biophysical components make their integration into global ecosystem simulation models very difficult18, where errors propagate and amplify uncertainty19. In this context, using expert knowledge and judgement with a rigorous and transparent methodology is one of the best paths toward progress20,21,22.

Fig. 1: Schematic representation of the sensitivity of ecosystem components to environmental change.
figure 1

Ecosystem components at critical thresholds (spheres located at the bottom of a slope) are exposed to the same persistent environmental change (indicated by the red arrows). The sensitivity of the components to change is illustrated by the time required (response time) to shift from an initial to a resulting state. Response times can be independent and homogenous (A1) or heterogenous (A2) among coexisting components. Interconnections between components (a sphere pulling other components) can reduce the heterogeneity in response times among coexisting components (A3). Response times can be homogenous across latitude (B1), homogenous between coexisting components but heterogenous across latitudes (B2), or heterogenous between components and across latitudes (B3). Note that definitions for specific terms are provided in the general glossary (Supplementary information).

Here, we provide a conceptual approach based on structured elicitation of expert judgement to quantify the relative sensitivity of the interacting components of northern ecosystems across a broad latitudinal gradient (a Glossary of key terms is provided in the Supplementary Notes). Expert judgement is an increasingly used methodology to tackle complex ecological questions at broad spatial scales23,24. Channelling the complementary expertise of specialists with decades-long first-hand experience of northern regions, we assess where and at what rate state shifts of ecosystem components are most likely to occur under warmer conditions along a ~100 km-wide, 3700 km latitudinal transect that extends from the boreal forest biome to the High Arctic (Fig. 2). This large study area located in eastern North America is among the last remaining wildernesses and the ecosystem states found along the entire latitudinal transect are primarily the result of natural postglacial processes, with limited direct anthropogenic impacts. While considering cascading effects and interactions between components, we assess the heterogeneity in response times under the scenario of a step change of a 5 °C increase in mean annual air temperature. This increase is well within the range of current warming scenarios for the Canadian Arctic (3–10 °C by the end of the 21st century relative to the 1986–2005 reference period; shared Socioeconomic Pathway SSP1-2.6 and SSP5-8.525).

Fig. 2: Location, extent and characteristics of the study region.
figure 2

The study region is a 100 km-wide strip of land covering a vast latitudinal gradient (3700 km) in northeastern North America. Mean annual air temperature along the transect currently ranges from −25–0 °C from north to south (see Supplementary Fig. 1 for July temperatures), while annual precipitation ranges from 100–200 to 800–1000 mm. Simplified illustrations of the current 28 contrasting states observed along the transect for six ecosystem components (permafrost, peatlands, lakes, winter snowpack, vegetation and endothermic vertebrate assemblages) are presented alongside the map. Detailed illustrations are available in Figs. 38. Black lines represent the southern limit (± 2.5 degrees of latitude) of each state identified along the latitudinal transect. The southern limit of the southernmost state for each component lies outside of the transect, and thus is not shown on this figure.

Results

The expert elicitation process identified a total of 28 contrasting states along the latitudinal transect (Fig. 2 and Table 1) for six key ecosystem components (permafrost, peatlands, lakes, winter snowpack, vegetation, and endothermic terrestrial vertebrates). For each ecosystem component, a minimum number of contrasting states (between 4 and 6) characterized by marked differences in structure and function were identified (implying that a state shift would strongly affect ecosystem properties). All states result from natural processes and are described in detail in Figs. 38 and Supplementary Note 1. The southernmost modern-day position of contrasting states along the transect (i.e., critical thresholds) were mostly located in the southern half of the transect (16 out of 22 positions are located in northern Québec; Fig. 9). We found that multiple critical thresholds emerged in the Subarctic region, between 55°N and 60°N. All ecosystem components were at critical thresholds within this region, characterized by the southern limit of continuous permafrost and of the shrub tundra. In the High Arctic, fewer thresholds emerged, and most were located slightly below 75°N, which corresponds to the southern limit of polar desert vegetation and intermittent ice cover lakes.

Table 1 Main characteristics used to discriminate contrasting states for six key ecosystem components of terrestrial landscapes along a ~100 km wide, 3700 km latitudinal transect across northeastern North America and examples of state shift impacts on ecosystem properties
Fig. 3: Permafrost contrasting states found along the transect.
figure 3

Permafrost is defined as rock or soil that remains below 0°C for at least two consecutive years. Along the latitudinal transect, it ranges from extensive continuous thick permafrost to relict permafrost located in isolated underground patches. Permafrost contrasting states are primarily based on the extent of frozen ground across the landscape (30 km X 30 km). States are also described using the characteristics (temperature, permafrost thickness, water and ice contents) of their different layers (active layer, perennially frozen ground and taliks – see Supplementary Note 1 for definitions).

Fig. 4: Peatland contrasting states found along the transect.
figure 4

Peatlands are terrestrial ecosystems in which waterlogged conditions prevent plant material from fully decomposing. Along the latitudinal transect, they range from polygonal peatlands shaped by the dynamics of ice-wedge networks to peatlands completely devoid of permafrost. The contrasting peatland states are defined by their surface morphology, the water and permafrost table dynamics, including the thickness of the active layer, as well as vegetation communities (see Supplementary Note 1 for definitions).

Fig. 5: Lake contrasting states found along the transect.
figure 5

Along the latitudinal transect, freshwater lakes (the typical lake defined for the study region has an average area of 24 ha and a mean depth of 4.0 m – see Supplementary Note 1) range from perennially frozen clear lakes to lakes with seasonal ice cover and high dissolved organic matter content. Lake contrasting states are defined by the duration of ice cover, stratification of the water column and mixing regime, oxygen and light availability, as well as the dissolved organic matter content (given as dissolved organic carbon, or DOC).

Fig. 6: Snowpack contrasting states found along the transect.
figure 6

Winter snowpack covers the study region for 6 to 10 months of the year. Snowpack contrasting states are described based on the depth of snow and snowpack stratigraphy, describing ice crystal type, grain size and hardness of the different snow layers of the snowpack in a landscape (30 |× 30 km) characterized by a flat terrain (see Supplementary Note 1 for definitions). Along the transect, the snowpack ranges from a two-layer poorly insulating, thin, hard and dense polar desert snowpack to a multi-layer, very thick and insulating boreal forest snowpack.

Fig. 7: Vegetation contrasting states found along the transect.
figure 7

Along the latitudinal transect, vegetation ranges from very sparse and low polar desert vegetation to dense and tall, closed canopy boreal forest vegetation. Vegetation contrasting states are defined by the dominant structural characteristics of vegetation cover, namely height and growth form (herbs, forbs, shrubs, coniferous and deciduous trees), found on a landscape (30 × 30 km) characterized by a flat terrain (see Supplementary Note 1 for definitions).

Fig. 8: Endothermic vertebrates contrasting states found along the transect.
figure 8

Along the latitudinal transect, endothermic terrestrial vertebrate communities range from an assemblage of few functional groups in the High-Arctic to a highly functionally diverse boreal vertebrate community. Contrasting states are defined by the occurrence and relative abundance of mammalian and avian functional groups present and active in the landscape (100 km x 100 km, a scale that encompasses summer home ranges of most taxa) throughout the annual cycle. Functional traits known to affect species interaction strength, dynamics of vertebrate communities and ecosystem processes were used to characterise functional groups (e.g., body size, nest placement, diet, hunting mode and sociality (see Supplementary Note 1 for definitions and details on functional groups and associated pictograms).

Fig. 9: Sensitivity of ecosystem components at critical thresholds.
figure 9

Expert assessment of the sensitivity to state shifts of ecosystem components at critical thresholds identified along a 3700 km latitudinal transect in northeastern North America. The coloured circles show contrasting states (initial and resulting) at each critical threshold. The response time is the estimated number of years (log scale) necessary for a component to shift from an initial to a resulting state under a warmer climate (i.e., a step change of 5 °C increase in annual air temperature). The southernmost modern-day ___location of a given contrasting state is its critical threshold. Overlaid spheres with letters indicate that a warming-induced change of a given component (L: lakes, P: permafrost, V: vegetation) is the key determinant of the response time of another component. Sensitivity estimates of all components are combined in the right panel, where locations of critical thresholds are vertically jittered for visualisation purposes.

The sensitivity of ecosystem components at critical thresholds, assessed using the time required to shift from an initial to a resulting state at the landscape scale (30 × 30 km, except for vertebrates at 100 × 100 km), ranged from high (1–10 years) to low ( >100 years). Several factors, such as precipitation, topography, soil characteristics, solar radiation, and wind exposure can modulate the response time of components to a step-change warming (i.e. increase or decrease response time within the assessed interval, see Table 2). The distribution of high and low sensitivities was not random among the observed critical threshold locations (Fig. 9). The mean latitude of critical thresholds associated with a low sensitivity was lower than randomly expected (p = 0.005), while the opposite pattern was found for critical thresholds associated with a high sensitivity (mean latitude higher than randomly expected (p = 0.002, Supplementary Figs. 24). Hence, even if low-latitude subarctic and boreal ecosystems are characterized by numerous critical thresholds, they should experience slower state shifts than some higher-latitude ecosystems close to critical thresholds and exposed to the same temperature increase.

Table 2 Key determinants of the sensitivity of the contrasting states of each ecosystem component to a hypothetical warming scenario assessed by expert knowledge

The non-random latitudinal distribution of high and low sensitivities likely reflects the predominance of cascading effects and strong interconnections among ecosystem components in the boreal and subarctic regions, as well as the relative effect of abiotic or biotic factors on response time along the latitudinal gradient (Fig. 9 and Table 2). The sensitivity of ecosystem components was often determined, at least partly, by the response time of another component (12 out of 22 state shifts). In most of these cases (8 out of 12), vegetation response time was a key determinant of the response time of other components such as permafrost, lakes, snowpack, and endothermic vertebrates. This underlines the ubiquitous cascading effects triggered by vegetation change, especially in the subarctic and boreal regions, where marked changes in plant cover and vertical structure occur. It also indicates that the relatively low sensitivity of a biological component can drive lower sensitivities in other ecosystem components and increase the coherence in the response times of coexisting components in some cases.

When vegetation was not a key determinant of the sensitivity of other components (9 out of 17 state shifts), response times were mainly determined by the rate of warming-induced changes in the cryosphere, as in the case of all peatland state shifts driven by permafrost degradation (Table 2). Changes in the cryosphere can sometimes trigger quick state shifts in ecosystem components, as shown by lakes at high latitudes, which were the most sensitive components of northern landscapes to warming in our study. For instance, a high arctic perennially ice-covered lake should shift to an intermittent ice cover lake over a relatively short timescale (1–10 years). On the other hand, when exposed to the same temperature increase, a seasonally ice-covered lake located at lower latitudes should shift to a contrasting state at a slower rate (10 to >100 years) because the response time is partly conditioned by changes in catchment vegetation that increase dissolved organic matter content (Fig. 9 and Table 2). Overall, the cascading effects of the most influential ecosystem components, vegetation and permafrost, appear to strongly drive the observed critical threshold locations and the heterogeneity in response times along the latitudinal gradient. By quantifying the effects of local environmental conditions that can modulate the response time of permafrost and vegetation to warming, especially precipitation, topography, and soil characteristics (Table 2), we could thus greatly enhance the accuracy of estimates regarding the sensitivity of northern ecosystems to climate warming.

Discussion

The identification and comprehensive understanding of regions ecologically sensitive to warming represent an urgent research priority on a global scale1,21. Our expert-based approach highlights clusters of critical thresholds and the strong heterogeneity in response times of ecosystem components to warming and associated cascading effects along a vast latitudinal transect of global significance. Our analysis identifies not only regions that are likely the most susceptible to state shifts, but also provides insights into which ones should most rapidly shift to contrasting states under warmer climates, and why. Theoretical studies suggest that some systems close to critical thresholds could be more sensitive to external perturbations17,26. However, our study reveals that terrestrial northern ecosystems characterized by multiple biophysical components at critical thresholds could be relatively resistant to future warming due to the lower sensitivities of key ecosystem components that strongly affect the response times of others.

The strong heterogeneity in the sensitivity to warming across latitudes and between coexisting components highlighted in our study should contribute to asynchronous expansion or contraction of the geographic ranges of different states of ecosystem components and translate to novel restructuring of ecosystems. For instance, in a High-Arctic landscape, the remarkably short response time of intermittent ice cover lakes may result in faster northward progression of seasonal ice-covered regimes compared to less sensitive coexisting component states, such as polar tundra vegetation, polar tundra winter snowpack, or polar tundra continuous permafrost. In some regions, the northward expansion of some ecosystem component states can be strongly influenced by human activities and by the presence of human-made or ecological barriers27,28. To better predict the northward progression of certain states, we would need to consider these effects in addition to the different warming scenarios and relative sensitivity of biophysical components.

A lack of coherence in the response times of biophysical ecosystem components to warming is expected to generate heterogeneous trends of biodiversity change across space and time28,29,30. The latitudinal southward decrease in sensitivity outlined in our study is consistent with the southward decrease in biodiversity response to climate change observed across northern ecosystems31. Moreover, the widespread northward (and alpine) advance of some terrestrial ecological boundaries appears much slower than the current warming rate, especially at lower latitudes32,33,34,35. Our results suggest that state shifts of certain terrestrial environments and their associated freshwater systems at lower latitudes are mainly driven by the response times of biological components with relatively low sensitivities to warming. However, in less productive ecosystems where there is limited or no erect woody vegetation, biophysical systems can transition swiftly to contrasting states if the response times of their components are primarily governed by physical processes, such as increasing ground and water temperatures.

Northern landscapes are trending away from states experienced during the 20th century and arguably from the conditions that prevailed over the past several millennia3,36. Under a warmer climate, it is possible that entirely novel states will emerge, and several climate variables could modulate the response times of ecosystem components. For example, large changes in precipitation (from snow- to rain-dominated; more frequent rain-on-snow events) are likely to accompany future warming trends37,38, amplifying the thermal effects on ecosystem change. We recognize that our results partly depend on the criteria used by experts to define biophysical contrasting states. As such, refining state definitions, adding more biophysical components and repeating the same analysis with different experts and climate change scenarios, and in various regions, would strengthen our conclusions.

Assessing and mapping sensitivity at large spatial scales is essential for predicting where and how ecosystem services will change in northern ecosystems, and for defining implications for indigenous communities who have been part of and dependent upon these ecosystems for millennia. Importantly, our study was not designed to accurately predict the future state of northern ecosystems, nor to determine exactly when state shifts will occur, but rather to assess the relative sensitivity and interplay of their different biophysical components. The multidisciplinary approach applied here, based on expert knowledge and judgement, represents a meaningful step forward in sensitivity assessment and allows the analysis of complex natural ecosystems across vast spatial scales. Our innovative approach could be applied to any region and provide the opportunity to map the sensitivity of multiple biophysical components across ecosystems of varying productivity and complexity, while also developing ways to refine response time estimates. Our study highlights the high degree of heterogeneity in response times between ecosystem components and a northward increase in sensitivity. The general lack of coherence in response times to warming between components and across latitudes, despite the presence of cascading effects, will likely impair the integrity of ecosystems by reshuffling the states of currently coexisting components and by generating heterogeneous rates of geographic range shifts. Such heterogeneity in sensitivity across ecosystems could ultimately lead to the restructuring, contraction, and expansion of biomes. In that context, the well-established concept of space for time substitution could lead to spurious results.

Methods

A supertransect across northern ecosystems

The study region is a ~100 km-wide band of land extending from the polar desert of northern Ellesmere and Ward Hunt (83°04’N), the northernmost islands of the Canadian Arctic Archipelago, to Matagami (49°45’N), located in the boreal forest biome of the James Bay region (Fig. 2). Multidisciplinary environmental studies have been conducted by numerous teams at several sites along this transect over the last 60 years39,40,41,42,43,44,45,46, providing a legacy of in-depth knowledge of the biotic and abiotic components that make up the landscapes of this vast northern region. The present condition of the various ecosystem components along the transect is partly the result of the late Pleistocene and Holocene history of the region47,48. Mean annual air temperature along the transect currently ranges from −25–0 °C, and annual precipitation from 100–200 mm to 800–1000 mm (according to CHELSA climate data, see Fig. 2)43,49. The study region encompasses other environmental variations with respect to hydrology, topography and geology. We restricted our analyses to low elevation areas (< 750 m), and excluded land currently covered by glaciers. Bedrock found along the entire transect belong to four broad geological provinces, mostly within the Canadian Shield (Superior and Churchill provinces) in the southern and central regions50. Hence, we assumed that geological differences across the study region have negligible effects on the ___location of critical thresholds (southern limits of contrasting states) and on the variation in sensitivity estimates. Moreover, it is a region with very sparse human population, negligible direct impact of human activities, and it is in one of the Earth’s last remaining wildernesses51,52. The current state of ecosystems found along the transect is thus essentially the result of natural postglacial processes. This makes the study region a global rarity that allowed us to focus on climatic impacts on the environment without the confounding effects of direct human pressures.

Expert elicitation

Expert judgement is often used to overcome uncertainty in complex scientific domains such as climate science22. The aim of our expert elicitation process was twofold: (1) to identify the ___location of critical latitudinal thresholds between ecosystem component states, and (2) to assess the sensitivity to warming of six key ecosystem components of northern landscapes: permafrost, peatlands, lakes, snowpack, vegetation and endothermic vertebrates. State shifts of these components can have strong impacts on ecosystem services (e.g., climate regulation, global biogeochemical cycles24,53) and on the culture, health and well-being of the indigenous communities inhabiting the study region (including food and water security, infrastructure, transportation, safety, health and traditional activities12,54,55,56).

We used structured expert elicitation inspired by the Delphi technique in which we combined multidisciplinary workshops with several rounds of expert assessments to ensure that all experts had the same interpretation of the questions and fully understood the aim of the project23,57. Our team is composed of 31 experts covering a wide range of disciplines (geomorphology, hydrology, nivology, limnology, plant and animal ecology and plant palaeoecology) and carrying out research across the entire study region. Experts were separated into 6 groups (one per component, 2 to 8 experts per group). Facilitators then addressed the following questions to each group: (1) What is the minimum number of contrasting states, characterized by marked differences in structure and function, that you can distinguish along the latitudinal transect (~100 km-wide transect)?, (2) At what latitude (± 2.5°) is the southern limit of each contrasting state located along the transect?, and (3) At the southern limit of a given state, how long would it take to shift to another, existing state (currently found along the transect) if the ecosystem was exposed to a step change of 5 °C increase in mean annual temperature? This last question was to be answered with a response time category reflecting High (1–10 years), Moderate (10–100 years), or Low (> 100 years) sensitivity. Facilitators also provided the following guidelines during the elicitation process: contrasting states and state changes had to be defined at the landscape scale (i.e. 30 km × 30 km, except for vertebrates at 100 × 100 km), state definitions had to be based on measurable characteristics with marked differences between states, and experts had to clearly describe the conditions for a state change as well as for the definition of a southern limit. Experts first identified contrasting states currently found along the transect. They provided detailed definitions for each contrasting state and justified their answers by providing detailed explanations (see Extended Data Figs. 16 for illustrations and Supplementary Note 1 for justifications and details). For each ecosystem component, experts used published and unpublished data to locate the southernmost ___location (± 2.5 latitudinal degrees) of all contrasting states along the latitudinal transect (e.g., the southern limit of continuous permafrost or the southern limit of shrub tundra vegetation; Fig. 2 and Supplementary Note 1). Ecosystems located at these southern limits were considered to be at critical thresholds, as ecosystems located further south (i.e., under warmer conditions) were characterized by a contrasting state.

Experts characterized the relative sensitivity of an ecosystem component at each critical threshold by assessing the order of magnitude of time necessary for a component located at the southern limit of its current state to shift to another existing contrasting state under a hypothetical warming scenario (3 categories of response time, in years: 1–10, 10–100, >100). Such broad time scales allowed experts to categorize response times of component state shifts to a marked temperature increase with greater confidence. Broad categories reflect the inherent variability associated with response time due, among other things, to variable local conditions and stochastic processes. We also used a log scale to reflect the increased uncertainty with increasing response time. Models predicting future warming scenarios indicate large annual air temperature changes at the end of this century for northernmost latitudes, with projected changes 2 to 4 times greater in Arctic regions relative to global trends4,25,58,59. In the Canadian Arctic, this corresponds to increases from 3 to 10°C (shared Socioeconomic Pathway SSP1-2.6 and SSP5-8.5) by the end of the 21st century (relative to the 1986-2005 reference period25). We thus used a step change of 5 °C increase in mean annual air temperature with respect to the 1986–2005 baseline period for the whole study area, and asked experts to categorize response times using time intervals on the log temporal scale (1–10 years, 10–100 years or >100 years). One of the key advantages of this approach is its capability to expose biophysical components to the same external forcing, enabling the comparison of relative sensitivity as a function of latitude and between coexisting components at the landscape scale. Experts considered known mechanisms, past changes and empirical studies, as well as feedbacks, cascading effects and interactions between ecosystem components. They also considered local conditions (e.g., topography, ground ice content, hydrology, precipitation) that can modulate response time. However, they did not fully integrate the impact of potential changes in multiple climate drivers. They nonetheless considered that temperature extremes and precipitation regimes would remain tied to the annual average temperature across the latitudinal gradient (i.e., a warmer climate is expected to increase annual precipitation). They also identified the main sources of uncertainty in response times that could be partly generated by synergistic links between some environmental variables with unknown trajectories that would ensue from warming (e.g., fire regime). Multiple response time categories were used in a few cases where the uncertainty was too broad.

When the initial round of questions to experts was completed, facilitators summarized the answers (number of contrasting states, definitions, southern limits, response times and justifications) and shared them among the groups. Five multidisciplinary workshops were organized for discussions between experts, who subsequently revised their answers and improved their justifications in smaller groups. Several rounds of revisions were made for each ecosystem component, during which experts were asked to comment, propose changes, and challenge the answers of other experts. A consensus on definitions, response time estimates, and justification was reached when facilitators obtained the approval of all experts (unanimous agreement among participants), independently. Our results thus reflect expert knowledge and judgement that is based on available model outputs, empirical data, and theoretical knowledge, but also unarticulated background information that enables the experts to transform their knowledge into a judgement22.

Statistical analysis

To test whether high, moderate and low sensitivities (response time intervals of 1–10, 10–100 and >100 years) were randomly distributed among the observed critical threshold locations (n = 22) we conducted permutation tests using the R software60. We first estimated the probability that the observed mean latitude of critical thresholds with low sensitivity (i.e., thresholds for which the maximum response time was >100 years, n = 8) was lower or higher than expected by chance. We computed the mean latitude for each possible permutation (n = 319,770). Using the resulting distribution, we calculated the p-value to determine if the observed mean latitude was unlikely to appear in a random situation (Supplementary Fig. 2). Following the same method, we estimated the probability that the observed mean latitude of critical thresholds with high and moderate sensitivity (i.e., critical thresholds for which the minimum response time was 1–10 years, n = 3; 1540 permutations, and critical thresholds with a response time of 10–100 years, n = 14; 319,770 permutations), was lower or higher than expected by chance (Supplementary Figs. 3,4).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.