Approaching a thermal tipping point in the Eurasian boreal forest at its southern margin

Climate change is increasing the intensity and frequency of extreme heat events. Ecological responses to extreme heat will depend on vegetation physiology and thermal tolerance. Here we report that Larix sibirica, a foundation species across boreal Eurasia, is vulnerable to extreme heat at its southern range margin due to its low thermal tolerance (T crit of photosynthesis: ~ 37 – 48 °C). Projections from CMIP6 Earth System Models (ESMs) suggest that leaf temperatures might exceed the 25 th percentile of Larix sibirica ’ s T crit by two to three days per year within the next two to three decades (by 2050) under high emission scenarios (SSP3-7.0 and SSP5-8.5). This degree of warming will threaten the biome ’ s continued ability to assimilate and sequester carbon. This work highlights that under high emission trajectories we may approach an abrupt ecological tipping point in southern boreal Eurasian forests substantially sooner than ESM estimates that do not consider plant thermal tolerance traits.

T he circumpolar boreal forest is the most extensive biome covering ~15 million km 2 and accounts for nearly a third of the 450Gt terrestrial carbon stock [1][2][3] .Boreal carbon resides in three major pools: peatlands, soils and forest biomass 4 .All pools are sensitive to rising temperatures from human-caused climate change due to direct impacts of warming and indirect impacts such as permafrost thaw, increased wildfire, and insect activity [3][4][5] .Two-thirds of the circumpolar boreal biome lies in Eurasia and North Mongolia, defined here as a ~1 million km 2 region between 95-114°E and 46-52°N, is located at the southern margin of the Eurasian boreal biome (Fig. 1a).The presence of boreal forest in North Mongolia can be attributed to its relatively cool mesic climate with mean annual precipitation generally exceeding 350 mm/year and mean annual temperatures between −5 and 0 °C (Fig. 1b and Supplementary Fig. 1a) 6,7 .Further south, the dominant ecosystem rapidly transitions into grasslandsteppe due to warmer and drier climate conditions 8 .
Boreal regions can host stable but competing forest, shrub, and grassland-steppe ecosystems as a result of various landatmosphere feedbacks such as albedo, fire, and insects 6,7,9,10 .Over the past four decades since the mid-1980s, accelerated warming and associated drying have decreased forest productivity and shifted the southern boreal forest margin northward, consistent with early stages of a projected ecosystem transition towards grass and shrublands 11,12 .Under continued global change, many Earth System Models (ESMs) predict the exceedance of a 'tipping point' in the southern boreal forest, a critical threshold at which a small perturbation in forcing (e.g., temperature) can fundamentally and irreversibly alter the ecosystem state at regional to sub-continental spatial scales (order of several hundred kilometres) over human timescales 13 .ESMs forecast synchronised southern margin boreal forest dieback to commence at ~1.5 C°of global warming, become widespread by ~3.5 °C, and exceed a biome-wide tipping point at ~4 °C of global warming relative to 1850 C.E pre-industrial levels 9 .Upon crossing this tipping point, southern boreal forest dieback is projected to occur over ~50-100 year timescales and result in up to ~52GtC of potential emissions 9,13 .Notably, North Mongolia and the circumpolar boreal forest have already warmed three times faster than global temperature due to Arctic amplification 10,[14][15][16][17][18] (Fig. 1c, d and Supplementary Fig. 2), with negative impacts on high-elevation Altai permafrost, tree growth 19 , and pastoral nomadic herding, the latter a traditional source of livelihood for a third of all Mongolians 8,[20][21][22] .ESMs participating in the sixth phase of the Coupled Model Intercomparison Project (CMIP6 23 ) suggest that warming will persist and intensify through the 21st century without compensating increases in precipitation (Fig. 1c, d and Supplementary Figs. 2, 3).in the context of North Mongolia (red rectangle), our study site (103.17°E,49.92°N, red star), and the circumpolar boreal biome (olive shading) 96,97 b Mean annual precipitation in Mongolia (1979-2021, CRU Ts 4.06 93 ).c Mean summer June-July-August (JJA) temperature for North Mongolia and the globe (area-weighted) between 1950-2021 93,94 , along with CMIP6 projected ensemble median 'historical' (1950-2014) and 'future' (2015-2100) temperature for both regions (22 models, 54 ensemble members, Supplementary Table 1) relative to their 1961-1990 mean (dashed zero line).Projections are derived using four shared socioeconomic pathways 95,98 (SSP1-2.6;SSP2-4.5;SSP3-7.0; and SSP5-8.5).The ensemble median weights models equally, accounting for variable ensemble member sizes across models.Further, it represents the common warming signal across models and is not expected to replicate the range of interannual variance of the observational data, unlike an individual ensemble member.Future projections are smoothed with a running 5-year mean (Supplementary Fig. 2 shows the model spread Rising mean temperatures are increasing the frequency and intensity of extreme heat events (i.e. the maximum temperature of the hottest day, T air-xx ) [26][27][28] .Changes in regional extreme temperature are expected to be greater than global mean temperature by a factor of 1.5 due to higher climate variability at local scales, differences in specific heat capacities of land and ocean, and feedbacks resulting from decreases in snow and soil moisture 27 .Further, while the intensity of warm extremes is expected to scale linearly with emissions and the degree of global warming, their frequency is expected to scale non-linearly 29 .The highest increase in T air-xx magnitude and frequency is projected to occur in mid-latitude regions, including North Mongolia 29 .These changes will have implications for forest ecosystems since all plant metabolic processes, including photosynthetic carbon assimilation, respiration and growth, are negatively impacted by temperature stress [30][31][32][33] .Photosystems I and II are two protein complexes that play a fundamental role in carbon assimilation and are located in the thylakoid membranes of plant leaves, where light energy is converted to chemical energy during photosynthesis 34 .The critical temperature of photosystem II (PSII) disruption (T crit ) represents a 'point of no return' at which leaves start to become irreversibly damaged and photosynthetic electron transport ceases 31,32,[35][36][37][38] .Sufficient accumulated damage to plants at regional scales can cause forest dieback and act as a trigger for abrupt ecosystem transition 32,35 .Extreme heat in particular has the potential to cause irreversible damage to plant photosynthetic apparatus since leaf temperature can exceed air temperature by +5-20°C depending on radiation load, transpiration, and wind speed [31][32][33]35,36,39 .
The importance of boreal forests in the global carbon cycle and the increasing probability of exposure to damaging high temperatures call for an urgent assessment of plant physiological thresholds to predict ecological change under a warming climate.However, most boreal thermal tolerance research has focussed on cold and not heat tolerance 29,32 despite boreal plants being particularly vulnerable to warming due to their low thermal safety margins 31 , defined as the difference between T crit and the maximal leaf temperature (T leaf-xx ) 40 .To the best of our knowledge, no characterisation of T crit exists for boreal Eurasia or across Asian forest biomes, more generally 32 .Additionally, most largescale ESMs currently do not include a warm-tolerance threshold in their representation of the temperature dependence of plant carbon assimilation 31 .Exceedance of this warm-tolerance threshold due to climate change can act as a trigger for abrupt ecosystem transition, substantially sooner than the ~50-100-year timescales over which the southern boreal forest transition is currently projected to occur 9,13 .
Here we apply a trait-based vulnerability assessment approach 41,42 to evaluate whether climate change may cause temperature extremes in excess of plant T crit at the southern margin of the Eurasian boreal forest.To do so we first comprehensively characterise the physiological responses of five dominant species in boreal North Mongolia by examining their chlorophyll fluorescence parameters, leaf traits, foliar nutrients, and foliar stable isotopes at a study site located in the Tarvagtai Valley, Bulgan in North Mongolia (red star in Fig. 1a, b).We then evaluate their thermal tolerance traits against CMIP6 ESMderived temperature projections through the end of the twentyfirst century.Species studied here include Siberian larch (Larix sibirica), silver birch (Betula platyphylla), Siberian elm (Ulmus pumila), Eurasian aspen (Populus tremula), and willow (Salix spp.) 43 (Supplementary Fig. 4).We place special emphasis on Larix sibirica since it is a foundation tree species across boreal Eurasia accounting for a third of the biome's total biomass [44][45][46] , and 80% of all Mongolian forest biomass where it provides important ecosystem services such as fuel wood and a productive grazing habitat for livestock 47 (Fig. 1a and Supplementary Fig. 5).

Results
Higher photosynthetic capacity but faster activation of photoprotection in Larix sibirica.Light energy reaching plant photosystems follows one of three main pathways 48 .Under favourable conditions, most light energy is used for photosynthetic electron transport (known as photochemistry), while excess energy is dissipated as heat (non-photochemical quenching, NPQ), and a small fraction (0-3%) is remitted as chlorophyll fluorescence 48 .We used Pulse Amplitude Modulated (PAM) chlorophyll fluorescence to estimate the photosynthetic capacity and performance of dark-adapted tree foliage 34,49 .We found that Larix sibirica converted a greater proportion of absorbed light energy to photochemistry as indicated by its higher maximum photochemical yield of PSII (F v /F m ) relative to all species (Supplementary Fig. 6a, p < 0.05 except Salix spp., pairwise nonparametric Wilcoxon signed-rank tests 50 ).
Rapid Light Curves (RLCs 48 ) indicated that Larix sibirica consistently achieved higher electron transport rates (ETR) for the same level of actinic irradiance (or photosynthetic photon flux density, PPFD) (Fig. 2a).This pattern of inter-species differences, with the highest values for Larix sibirica, was also observed for quantum efficiency (α), maximum ETR (ETR max ), and minimum saturating irradiance (I k ) (Supplementary Fig. 6b-d).There were no consistent differences in ETR for sun and shade leaves (Supplementary Fig. 7).RLCs for photochemistry expressed by the effective quantum yield of PSII (Φ PSII ) and for heat dissipation (NPQ) showed that as irradiance (PAR, photosynthetically active radiation) increased, Φ PSII decreased exponentially while NPQ increased following a saturation curve for all five species 51 (Supplementary Fig. 8).
All species experienced a dynamic reallocation of absorbed light energy from photochemistry to NPQ as irradiance intensity increased, as shown by the negative relationships between Φ PSII and NPQ (Fig. 2b and Supplementary Fig. 9).However, Larix sibirica (and Ulmus pumila) exhibited stronger and more negative linear relationships between Φ PSII and NPQ indicating higher vulnerability to damage at increasing light levels and a faster activation of photoprotection via NPQ.Φ PSII is a measure of the proportion of light energy absorbed by PSII used in photochemistry, while NPQ is indicative of the amount of absorbed light energy dissipated by changes in the xanthophyll cycle deepoxidation state (i.e., not used for photochemistry or re-emitted as fluorescence) 48,51 .This dynamic reallocation is associated with increased photoprotection to prevent damage to the photosynthetic apparatus at high light 52 .Nevertheless, Larix sibirica maintained a higher ETR and Φ PSII and a lower NPQ relative to all species at all irradiance levels (Supplementary Figs. 8, 9).
Larix sibiricaʼs higher abundance and biomass across the landscape 44,46 may therefore partially be explained by how it outperforms co-located species in carbon assimilation under current climate conditions.An ETR of 4 µmol/m 2 s (y-axis in Fig. 2a) is theoretically equivalent to a gross photosynthesis rate of 1 µmol/m 2 s (ref. 51).However, while ETR and CO 2 assimilation are generally proportional, their relationship can be non-linear for C 3 plants 51 .CO 2 assimilation is also a function of the carboxylative activity of the enzyme RuBisCO, which is in-turn influenced by stomatal and mesophyll conductance 34 .Additionally, the oxygenation by RuBisCO instead of carboxylation (photorespiration), cyclic electron transport, and nitrate reduction can be alternative electron sinks to carbon assimilation 34 .Consequently, we analysed the foliar chemistry of all species to further examine Larix sibirica's higher photosynthetic capacity 53 .
Enriched stable carbon isotopic ratios in Larix sibirica as a result of high photosynthetic rates.Foliar stable carbon (δ 13 C) and oxygen (δ 18 O) isotopic ratios showed a consistent pattern of being most enriched for Larix sibirica and most depleted for Populus tremula (δ 13 C) and Salix spp.(δ 18 O) (Fig. 2c and Supplementary Fig. 10a).However, foliar δ 13 C and δ 18 O were not correlated in Larix sibirica (Spearman r = −0.06,Fig. 2d).Foliar δ 13 C is influenced by the rate of CO 2 fixation by the enzyme RuBisCo (photosynthetic rate, A) and stomatal conductance (g s ), while foliar δ 18 O is primarily determined by soil source water and evaporative enrichment in the leaf during transpiration 54,55 .Since stomatal conductance (g s ) can be a determinant of both δ 13 C and δ 18 O, exploring δ 18 O in conjunction with δ 13 C (i.e., the dualisotope method) can shed some light on whether changes in δ 13 C are related with the rate of photosynthesis (A) or stomatal conductance (g s ) 56 .Consequently, as leaf δ 18 O is not influenced by the rate of photosynthesis (A), we interpret that Larix sibirica's enriched leaf δ 18 O reflects site-level microclimate differences (Supplementary Fig. 4).Despite being a mesic region overall, δ 18 O was positively correlated with elevation (Spearman r = 0.39-0.42,p < 0.01) in correspondence with a gradient of better-drained montane soils that undergo higher evaporative enrichment than river valleys (Supplementary Fig. 10b).The lack of relationship between foliar δ 13 C and δ 18 O, together with likely microclimatic causes for the δ 18 O enrichment, suggests that the enriched δ 13 C in Larix sibirica is likely related to higher rates of carbon assimilation 56,57 .
Photosynthetic capacity can also sometimes be associated with higher nitrogen availability since nitrogen is a building block for the photosynthetic enzyme RuBisCo 7,58,59 .However, similar foliar δ 15 N across all species suggests a lack of site or specieslevel differences in nitrogen availability, source (nitrate|NO − 3 cf.ammonia|NH + 4 ), or root mycorrhizal associations 57,60,61 (Supplementary Fig. 10c).Supporting this result, Larix sibirica leaves had lower foliar nitrogen [N%] and higher carbon [C%] concentrations resulting in higher carbon to nitrogen (C:N) ratios (Supplementary Fig. 11).
Therefore, the high photosynthetic performance in Larix sibirica is supported by several lines of evidence such as higher ETRs, quantum efficiency (α), F v /F m , and enriched foliar δ 13 C, but does not seem related to higher nitrogen availability 46 .This higher performance of Larix sibrica in relation to other species is likely a species-specific trait that facilitates its current role as a foundation species in the North Mongolian and Eurasian boreal forest.
Larix sibirica's low T crit may make it more vulnerable to damage from high temperatures.Larix sibirica T crit is 3.5-6.7 °C lower than other North Mongolian species (Fig. 3 and Supplementary Figs. 12, 13).While warm temperatures can decrease plant productivity, warm extremes play a more important role in determining plant mortality, survival, photosynthetic performance, and even plant evolution 30,32 .At T crit , minimal chlorophyll-a fluorescence of PSII (F o ) begins to rise rapidly, indicating damage to thylakoid membrane where the PSII is located 37,62,63 (Supplementary Fig. 12).An evaluation of T crit of Larix sibrica in terms of observed and projected T air-xx under different SSPs suggests overlap between the lower threshold of the Fig. 2 Plant photosynthetic performance and foliar isotopes.a Electron transport rate (ETR, µmol-electrons/(m 2 s)) against photosynthetic photon flux density (PPFD, µmol photons/m 2 s).PPFD is the intensity of photosynthetically active radiation (PAR) in the ~400-700 nm spectral range that plants can use for photosynthesis 99 (full sunlight is generally between ~900-1500 μmol photons/m 2 s).See Supplementary Fig. 6 for quantum efficiency (α), maximum ETR (ETR max ), and minimum saturating irradiance (I k ).Error bars describe ± one standard error.The number of sampled leaves is noted in the legend.b.An inverse relationship between the effective yield of Photosystem II (Φ PSII ) and non-photochemical quenching (NPQ) with increasing PAR in Larix sibirica.c Foliar δ 13 C for all species showing that Larix sibirica is significantly enriched in δ 13 C (p < 0.05, Wilcoxon signed-rank test 50 ).Each box describes the 25th, 50th (i.e.median) and 75th percentiles of data and whiskers describe the range.Statistically significant differences are indicated by the lack of a shared letter label between species, while numbers describe the species mean.d Scatterplot between foliar δ 13 C and δ 18 O for all species, along with a loess fit and Spearman rank correlations.
current thermal tolerance range of Larix sibirica and T air-xx under SSP3-7.0 and SSP5-8.5 in the second half of this century (2050-2100) (Fig. 3a, b).These results indicate that Larix sibirica is more vulnerable to damage to its photosynthetic apparatus than other co-located species in North Mongolia and would be more frequently exposed to damaging summer temperatures exceeding its T crit as summer temperatures continue to warm through this century 32,37,53,63,64 (Fig. 1c and Supplementary Fig. 2a).
The T crit of a species may, however, acclimate with warming 31,32,65 .Acclimation in some species can occur seasonally, with lower T crit during cool winter months (in evergreen species) and higher T crit during summer months, or under exposure to warmer growing conditions 63,66 .In a global survey of temperature responses, T crit increased on average by 0.38 °C per °C increase in the mean maximum temperature of the warmest month 31 .This global increase in T crit with temperature is greater than the observed regional and seasonal acclimatisation of T crit 38,63,67 .We, therefore, applied a 0.38 °C increase in T crit per °C increase in the ensemble median T air-xx as an optimistic upperrange estimate of thermal acclimation of T crit for Larix sibirica (Supplementary Fig. 13b).Nevertheless, we note that the mechanisms of thermal acclimation are poorly understood 63,68 .It is uncertain whether T crit acclimates in all species, particularly in deciduous species such as those studied here, and whether there are limits to the degree of acclimation 63,69 .Additionally, the rate of increase in temperature due to climate change is projected to be greater than the rate of T crit acclimation (see the following section).Therefore, while thermal acclimation may potentially moderate some of the impacts of extreme T air-xx 68 it is likely that Larix sibirica may be more vulnerable to rising temperature if its T crit acclimates at a slower pace or negligible pace than assumed here.
Leaf temperature exceeds Larix sibirica's T crit under high emission scenarios by 2050.Leaf temperature (T leaf ) is a more relevant metric to evaluate plant thermal tolerance since T leaf can exceed air temperature by +5-20 °C31,32,38,40,70 .On diurnal timescales, photosynthesis and stomatal conductance generally decrease during the early afternoon due to high temperature and a high leaf-to-air vapour pressure gradient [71][72][73] .At longer daily to seasonal timescales, stomatal closure can occur during drought conditions 30 .Partial to full stomatal closure during the 'mid-day depression' of photosynthesis and during drought reduce plant transpiration and associated latent heat flux cooling and enhances leaf temperature above air temperature 30,72,74 .
Higher T leaf can be expected regionally in coming decades as a result of increased stomatal closure from enhanced VPD, T air , atmospheric CO 2 and the lack of projected increases in precipitation 24,74 (Figs. 3, 1c and Supplementary Figs.1-3).To overcome the lack of field measurements of T leaf we used ERA5 surface skin temperature (T skin ) as a proxy for T leaf (by setting T leaf equal to T skin ).Doing so allows us to estimate the leaf temperature (T leaf ) in relation to T air .T skin is the approximate theoretical temperature required to satisfy the surface energy balance 75 , and represents how hot an object would feel to touch since sun-exposed objects can become considerably warmer than surrounding 2 m air temperature 76 .The climatology of our site shows that T air , T skin, and precipitation peak during the summer JJA months (Supplementary Fig. 14).Further, T skin-xx is generally 1-6 °C warmer than T air-xx (Supplementary Fig. 15) with a ~4 °C difference between median T skin-xx and T air-xx .We apply this 4 °C offset to T air-xx projections as an estimate of T leaf-xx .This offset is conservative since it is calculated over a coarse spatiotemporal field (a 0.2°× 0.2°spatial resolution hourly reanalysis data), assumes an infinitesimally thin surface layer without thermal memory, and is considerably lower than field estimates 31,32,36 .Despite this, we find marked increases in the likelihood that T leaf-xx exceeds Larix sibirica's warming-acclimated T crit by 2050 under the higher emission scenarios SSP5-8.5 and SSP3-7.0,but not under lower emission scenarios SSP2-4.5 and SSP1-2.6 (Fig. 3c, d and Supplementary Figs. 16, 17).
Under the high-emission SSP5-8.5 and SSP3-7.0scenarios, starting around 2050, annual T leaf exceeds the 10th to 25th percentile of Larix sibirica T crit at least one day per year in 5-20% of CMIP6 simulations (Fig. 4a, b and Supplementary Fig. 18) and Larix sibirica's median T crit in 2-10% of the 56 CMIP6 ensemble members (Supplementary Fig. 19).Additionally, the 25th percentile of Larix sibirica's T crit is exceeded a median of 2 to 3 days per year in SSP5-8.5 and up to 2 days per year under SSP3-7.0around 2050 (Fig. 4c, d and Supplementary Fig. 19).Similar results are found for SSP5-8.5 and SSP3-7.0 using only one ensemble member per model and avoiding overweighting models with more than one ensemble member at the expense of lower sample size of model runs (Supplementary Fig. 20).
In addition to climate change-driven warming, annual T leaf and T air are also controlled by natural year-to-year climate variability.Therefore, although the exceedance T crit by T leaf in any given year should theoretically be a low-probability event, we show that exceedances of Larix sibirica's 25th percentile of T crit occur across 22 of 23 ESMs used (Supplementary Fig. 21).This suggests that there is potential for damage to Larix sibirica's photosynthetic apparatus in the next three decades by ~2050 under SSP5-8.5 and SSP3-7.0 and that these results are not biased by choice of CMIP6 models used in the presented analysis.However, we note that damage is unlikely under lower emission scenarios (SSP2-4.5 and SSP1-2.6;Supplementary Figs. 16, 17), and projected damage to Larix sibirica's photosynthetic apparatus generally remains less than 60% under the two high emission scenarios SSP5-8.5 and SSP3-7.0(Supplementary Fig. 18).

Discussion
Boreal ecosystems contain about one-third of the total terrestrial carbon stock and account for approximately one-third of the annual terrestrial carbon sink 2,77 .A decreased boreal forest carbon sink and release of current carbon stocks following large scale forest mortality can potentially contribute to a 'positive climate feedback' and exacerbate ongoing warming 78 , although we note that net radiative forcing will also depend on surface albedo change 9 .Further, boreal regions in the northern hemisphere have warmed considerably faster than the global average 29,79 .JJA temperatures in North Mongolia have warmed approximately three times the global average (1.81 °C cf.0.63 °C) over the past few decades (2001-2021 relative to 1961-1990) and are expected to continue to warm more rapidly than the globe under a range of low to high emission scenarios (Fig. 1 and Supplementary Fig. 2).
Extreme heat can cause permanent damage to a plant's photosynthetic apparatus and can reduce carbon assimilation by forests 31,32 .Damage to the photosynthetic apparatus can occur at the time scale of seconds to minutes following heat exposure 35 .Additionally, forests are highly sensitive to hot temperature extremes at the warm margins of a biome or species distribution range 41,42 .By applying a trait-based vulnerability assessment approach 41,42 , we found that Larix sibirica had a significantly lower T crit among other co-located species in the southern margin of the Eurasian boreal forest (Fig. 3b).Further, we project increases in the frequency, magnitude, and duration of exceedances of Larix sibirica's T crit using a suite of 23 CMIP6 ESMs (56 ensemble members) if greenhouse gas emissions are not reduced rapidly within the next two decades.T crit exceedances are projected to occur almost annually starting ~2050 under higher emission scenarios (SSP5-8.5 and SSP3-7.0),but not under lower emission scenarios (SSP2-4.5 and SSP1-2.6)(Figs. 3, 4).The rise in frequency and duration of heat extremes under high emission scenarios, measured by the number of days per year that T leaf exceeds T crit , is likely to be particularly deleterious to plant health.While trees can recover from a single stress event if they occur infrequently, frequent and repeated exposure to stress may lead to accumulated damage and exhaustion of carbohydrate reserves 39,80 .This can weaken plant resilience enough that even moderate stress might cause whole plant mortality 39,80 .Additionally, we project increasing aridity in the region due to a higher VPD caused by increasing temperature and a lack of a compensating increase in precipitation 24 .While we assume a conservative 4 °C offset between T air-xx and T leaf-xx we note that plants may be more vulnerable to heat damage in the future as stomatal closure due to enhanced drought stress and elevated atmospheric CO 2 can elevate leaf temperature 25,81 .
The T crit of a species has been hypothesised to be related to the temperature of the environment it evolved or is currently distributed in (evolutionary influence on T crit ), the temperature regime of the local site (phenotypic influence on T crit ) 32,63,65 , and soil moisture availability 81 .The differences in T crit observed here among the five species cannot however be explained by elevation differences alone (Supplementary Fig. 13) or micro-site differences in soil moisture.Among the four species for which distribution maps are available, Larix sibirica is the only species that is located at the southernmost margin of its distribution range at our study site 44,46 (Fig. 1a and Supplementary Fig. 4b).The remaining three species (Betula platyphylla, Ulmus pumila, and Populus tremula) are also distributed further southeast of our study site in warmer and drier localities (Supplementary Fig. 22) 82,83 .Thermal tolerance has also been shown to be higher for species in more humid sites since wetter soils may allow leaves to keep their stomata open for longer and benefit more from evaporative cooling at the leaf surface 81 .In this regard, while Larix sibirica's lower T crit relative to other species is unexpected at the site-scale based on its preference for drier micro-habitat conditions in better-drained soils away from the river valley (Supplementary Fig. 4a), it could be related to its more northerly distribution in mesic boreal Eurasia.
The lower T crit of Larix sibirica and Ulmus pumila relative to other species (Fig. 3b) is also supported by the stronger negative linear relationship between Φ PSII and NPQ observed for these two species (Supplementary Fig. 9).The relative rate of increase in NPQ compared to Φ PSII at moderate to high PAR values is greater for Larix sibirica and Ulmus pumila than other species which suggests the need for enhanced photoprotective mechanisms at lower light levels.This suggests that they are more 'vulnerable' to damage at increasing light.Other species generally show a nonlinear response between Φ PSII and NPQ and can maintain relatively low NPQ under moderate light levels.Nonetheless, Larix sibirica likely outperforms all species in photosynthetic performance as suggested by chlorophyll fluorescence and foliar stable carbon isotopic data (Fig. 2).This higher photosynthetic performance is not related to differences in nitrogen availability (Supplementary Figs. 10, 11) and is also likely a species-specific trait 46 that enables Larix sibirica to be highly successful in the Eurasian boreal forest under current climate conditions.
Since Larix sibirica is a foundation species across 30% of the Eurasian boreal forest and is the most common tree in Mongolia 47,84 , an increase in the frequency and intensity of extreme heat events in excess of its T crit may represent a thermal tolerance tipping point at the southern margin of the Eurasian boreal forest.Without rapid emission reductions within the next two to three decades prior to ~2050, our results suggest that the viability of Larix sibirica in the boreal forests of Mongolia will likely be threatened by repeated incidences of damage to its photosynthetic apparatus.Combined with the increased likelihood of additional disturbance agents such as drought, fire, insects, and pathogens 85 , such damage has the potential to result in large-scale mortality of Larix sibirica-dominated forests, as has already been observed in some locations (Supplementary Fig. 23), and a possible conversion of the ecosystem into a grassland-steppe biome 86 .These changes will likely have significant implications for the region's forestry and herding sectors.However, we do note that our projections of future T leaf-xx generally do not exceed the 60th percentile of Larix sibirica's T crit and never exceed its maximum T crit under either high emission scenario SSP5-8.5 and SSP3-7.0.This suggests a scenario where Larix sibirica individuals growing in favourable microclimate conditions (e.g., on cooler north-facing high-elevation slopes) might be able to escape substantial damage from the increasing incidences of extreme heat events.
Additionally, we also acknowledge a few limitations inherent in the sampling design of our study that may contribute to additional uncertainty.We only targeted mature dominant trees at one site in Northern Mongolia in one year, during the warmest part of the year (~5 species × 5 individuals × 2 samples/individual, in July-August 2019).Our results, therefore, may not adequately characterise the within-species variation of the measured traits (T crit , RLCs, foliar chemistry) both at the site and across boreal Eurasia, their dependence on stand age and density, their variability across years and during the course of the year (note that all species are deciduous).In this regard, we were primarily limited by the timeintensive nature of running over 50 RLCs and temperaturefluorescence curves that took more than one hour per sample in a very remote location without access to grid electricity.We conducted our sampling during the summer months of July and August since plant physiological parameters, including photosynthetic performance and warm thermal tolerance, exhibit seasonal patterns in temperate and boreal climates and generally reach their peak values during the summer months 69 .Tracking of plant physiological parameters through the active season will therefore aid in more fully evaluating plant stress tolerance, for example, to spring season heatwaves when plant T crit are at their lowest relative values to temperature, but are beyond the scope of this work.Additional plant ecophysiological analyses are therefore needed to more fully characterise the range of spatial and temporal plant thermal tolerance and physiological variability across boreal Eurasia to comprehensively assess the biome's vulnerability to climate change.Further, in our study, we focus on the frequencies of extreme heat events (T leaf in excess of T crit ) and not on the absolute magnitude of future T leaf itself.Downscaled and bias-corrected CMIP6 ESM output (e.g.ref. 87 ) can be used to better characterise changes in the magnitude of future extreme heat events (in °C) in the region, particularly in relation to the exceedance of plant T crit , though we note that they should not adjust the frequency of extreme heat events.

Conclusions
Analysis of tipping points in the southern boreal biome margin using ESMs have suggested that large-scale forest dieback will likely occur under high emission scenarios (SSP3-7.0 and SSP5-8.5)by the end of this century over timescales of ~50-100 years 9,13 .However, the relationships between forest mortality and disturbance agents such as fire, drought and insects are poorly characterised in ESMs 88 .Additionally, most ESMs currently represent the temperature dependence of plant carbon assimilation using Arrhenius-type functions that assume a thermal optimum for photosynthesis and a fully reversible decline in photosynthesis at high leaf temperature 35,89 , while in-situ studies show that even a few minutes of exposure of leaves to temperatures in excess of T crit can cause long-term damage to plant photosynthetic apparatus and result in tissue necrosis 35 .Consequently, it is possible that we may currently be underestimating the timing and rate of southern boreal biome ecosystem change 88 .
Here we identify potential mechanisms for an abrupt thermal tipping point in southern boreal Eurasia due to projected increases in the frequency and intensity of extreme heat events (T air-xx ) under continued climate change that may cause T leaf to exceed T crit .Such exceedances, particularly if they occur in close succession 80 , may act as a trigger for regional-scale forest mortality on timescales of days to years.Though we note that eventual tree death will likely be caused by a combination of multiple stress factors such as hydraulic failure, depletion of carbon reserves, drought, insects, and pathogens, in addition to direct heat damage.We find that this tipping point may be exceeded as soon as 2050 under high-emission trajectories (SSP3-7.0 and SSP5-8.5) 90, significantly sooner than prior ESM estimates that predict ecosystem transitions to occur at the end of this century over decadal to centennial timescales.Incorporating plant T crit in the next generation of 'trait-based' vegetation models and ESMs will therefore be important to accurately simulate the future terrestrial carbon cycle both globally and in the Eurasian boreal forest.

Methods
Study site and foliar sampling.Our study site is located at 103.17°E, 49.92°N in the Tarvagatai River valley in the Bulgan aimag (province) of north-central Mongolia and east of the Teshig soum.The Taravagatai Valley is located in the forest-steppe region of the Baikal Rift Zone and has a basal elevation of ~900 m above sea level 91 .The river Tarvagatai originates in the Khantai mountains and flows through this valley, eventually joining the Eg river and forming one of the largest rivers in Northern Mongolia 43 .The regional climate is strongly seasonal, with long and cold winters and short but warm summers (Supplementary Fig. 14) 43 .In total, 53 leaves were collected from a total of 25 trees (5 trees × 5 species).On average, we collected two leaves per tree (25 trees × 2 leaves/tree), with a sun leaf collected on the south-facing aspect and a shade leaf collected on a northfacing aspect.We collected leaves at multiple elevations within ~5-7 km of the study site on dominant tree individuals.
Rapid light curves (RLC), foliar chemistry and critical temperature of PSII disruption (T crit ).We made RLC measurements using the Junior PAM portable chlorophyll fluorometer (Heinz Walz GmbH).Each RLC consists of the fluorescence responses to nine different increasing actinic irradiances (i.e.light energy that can trigger plant photochemistry) ranging from 0 through 845 µmol photons/ m 2 s.The Junior PAM fluorometer is fitted with a 1.5 mm diameter fibre optic and a blue diode (485 ± 40 nm) with an attachable magnetic clamp on one end that holds down the leaf specimen to ~1 mm from the end of the optic fibre 51 .We measured PAM parameters such as the Electron Transport Rate (ETR), the quantum efficiency of photosynthesis (α), non-photochemical quenching (NPQ), etc., using the WINCONTROL-3 software.The apparent rate of photosynthetic electron transport of PSII (ETR) was obtained as: ETR ¼ Φ PSII Ã PPFD Ã 0:5 Ã 0:84 where, PPFD represents the photosynthetic photon flux density (in µmol photons/ m 2 s) of the applied actinic irradiance.The effective quantum yield of photosystem II (Φ PSII ) is calculated as (F m '-F)/F m ', where F m ' is the maximal fluorescence after each incremental light pulse and F is the fluorescence during light treatment.The scaling factor of 0.5 assumes equal excitation of both PSII and PSI, while the absorptance factor of 0.84 is an estimate of the fraction of incident light that is actually absorbed by PSII 48,92 .Prior to performing the RLCs, all leaves were wrapped in moist paper towels and were placed in the dark for a minimum of 2 h.This dark adaptation minimises the impact of reversible photoinhibition by permitting the dissipation of the electrochemical gradient across the thylakoid membrane 48 .However, to further drain electrons from the acceptor side of PSII and fully oxidise it we applied a 5.5 s far-red-light (6 µmol m −2 s −1 ) treatment prior to all measurements 48 .RLCs were conducted with leaf temperatures around 20 °C.Following an RLC, each leaf was allowed to rest in the dark under the magnetic leaf clip for ~5 min until minimal chlorophyll-a fluorescence levels stabilised, after which another 5.5 s far red-light treatment were applied on the leaves 53 .Finally, each leaf was heated from 20 °C to 60 °C at a rate of ~1 °C/min and minimal chlorophyll-a fluorescence was measured continuously using the Junior PAM device to produce fluorescence-temperature (FT) curves (Supplementary Fig. 12).T crit was determined by fitting a piecewise best-fit linear regression to each FT curve 31 .
We then dried each leaf for a minimum of 48 h at 60 °C in a drying oven and ground them in separate vials using a ball mill grinder.Following this, we determined leaf percent carbon (%C), percent nitrogen (%N), percent oxygen (%C) using an elemental analyser (ECS 4010, Costech Analytical Technologies, Inc., Valencia, California, USA), and leaf δ 13 C, δ 18 O and δ 18 O using an isotope ratio mass spectrometer (Delta PlusXP, Thermo Finnigan, Bremen, Germany).Climate and Earth system models.We derived climate data from three datasets, (i) CRU Ts 4.06 monthly data 93 , (ii) HadCRU5 94 and (iii) ERA5 75 .We downloaded the earth system model (ESM) output from the Coupled Model Intercomparison Project (CMIP6) archive at PANGEO (https://pangeo-data.github.io).CMIP6 includes a suite of model experiments organised and contributed from multiple modelling centres in support of the Sixth Assessment Report (AR6) of the Intergovernmental Panel on Climate Change (IPCC) 23,95 .Our analysis is restricted to models that had continuous ensemble members spanning the 'historical' (1950-2014 C.E.) and 'future' (2015-2100 C.E.) time periods to allow for a comparison between model simulated 'historical' and 'end-of-the-century' climate on a model-by-model basis.Further, we only used ESMs that had both daily and monthly model outputs available on PANGEO.This resulted in a subset of 23 models and 56 ensemble members (Supplementary Table 1).However, for analysis of monthly temperature data, only 22 models and 54 ensemble members are considered as the historical period temperature data was not available for MPI-ESM1-2-HR (footnote in Supplementary Table 1).
To calculate projections of summer JJA temperature and precipitation (Fig. 1 and Supplementary Figs. 2, 3) we first downloaded monthly 2 m near-surface air temperature (variable tas) and monthly precipitation flux (variable pr) for the globe and for North Mongolia (46-52°N, 95-114°E) from the CMIP6 archive.Next, we converted the monthly precipitation flux (kg/(m 2 .s))into monthly precipitation (in mm).We then calculated all climate series (ERA5 temperature, CRU Ts 4.06 precipitation, all CMIP6 ensemble members) for the globe and for North Mongolia relative to their 1961-1990 mean.Following this, we computed T air-xx as the maximum temperature of the hottest day of each year using ERA5 and CMIP6 data.To achieve this for ERA5, we first calculated the daily maximum T air using hourly ERA5 data, while for CMIP6 data, we directly downloaded daily maximum temperature data for each model ensemble member (variable tasmax).CMIP6 tasmax data were averaged across a 1.5 × 1.5-degree latitude-longitude grid box around our site for projections of T air-xx .To develop projections of T air-xx we first adjusted the mean T air-xx of each ensemble member between 1979 and 2014 to match the mean ERA5 T air-xx for the same period.We applied this 'mean correction' since 45 of the 56 CMIP6 ensemble members exhibited a 'warm bias' relative to ERA5 data (Supplementary Fig. 24).On average, CMIP6 projections of T air-xx would be 1.67 °C warmer without this mean value adjustment.

Fig. 1
Fig.1Climate trends and projections in North Mongolia and the globe.a Distribution of Larix sibirica (black hatching)44,46 in the context of North Mongolia (red rectangle), our study site (103.17°E,49.92°N, red star), and the circumpolar boreal biome (olive shading)96,97 b Mean annual precipitation in Mongolia (1979-2021, CRU Ts 4.0693 ).c Mean summer June-July-August (JJA) temperature for North Mongolia and the globe (area-weighted) between 1950-202193,94 , along with CMIP6 projected ensemble median 'historical' (1950-2014) and 'future' (2015-2100) temperature for both regions (22 models, 54 ensemble members, Supplementary Table1) relative to their 1961-1990 mean (dashed zero line).Projections are derived using four shared socioeconomic pathways95,98 (SSP1-2.6;SSP2-4.5;SSP3-7.0; and SSP5-8.5).The ensemble median weights models equally, accounting for variable ensemble member sizes across models.Further, it represents the common warming signal across models and is not expected to replicate the range of interannual variance of the observational data, unlike an individual ensemble member.Future projections are smoothed with a running 5-year mean (Supplementary Fig.2shows the model spread).d Observed warming in North Mongolia and the globe between 2001-2021 (OBS.) and at the end-of-thecentury (2081-2100) using four CMIP6 SSPs (relative to 1961-1990).CMIP6 boxplots represent spread across 54 ensemble members.The red dot and the text at the bottom describe the median warming for observational data and ensemble median warming for CMIP6.Observed JJA warming in North Mongolia is nearly thrice the global average (1.82 °C cf.0.63 °C, 2001-2021 cf.1961-1990).The region is projected to warm ~1.7 times the global average through the 21 st century, regardless of the SSP used.Same y-axis label in c and d with different ranges.Each box describes the 25th, 50th (i.e., median) and 75th percentiles of data and whiskers describe the range.

Fig. 3
Fig. 3 Maximum temperature of the hottest day (T air-xx ) and hottest leaf temperature (T leaf-xx ) in relation to the critical temperature of photosystem II disruption (T crit ) for all five study species.a ERA5 T air-xx between 1959-2021 (in blue, averaged across a 0.2 × 0.2-degree grid box around our study site: 103.17°E, 49.92°N) along with CMIP6 projected ensemble median 'historical' (1950-2014) and 'future' (2015-2100) T air-xx using 4 four shared socioeconomic pathways (SSPs) (23 models, 56 ensemble members, averaged across a 1.5 × 1.5-degree grid box around our site).Horizontal and vertical lines represent T air-xx of 40 °C and the year 2050, respectively.b Range of variability in T air-xx between 2050-2100 under different SSPs (red dot: ensemble median) and measurements of the T crit of all five study species measured in the field in 2019.Ensemble median T air-xx and median T crit are described below each boxplot.The lack of a shared common letter label between Larix sibirica and all other North Mongolian species indicates that it has a significantly lower T crit (p < 0.05 using a non-parametric Wilcoxon signed-rank test 50 ).Larix sibirica's measured thermal tolerance range [T crit-min , T crit-max ] is 37.15-47.96°C (mean 43.25 °C, median 43.18 °C).Each box describes the 25th, 50th (i.e., median) and 75th percentiles of data and whiskers describe the range.c, d Projections of future T leaf-xx under SSP5-8.5 (c) and SSP3-7.0(d) in relation to projections of Larix sibirica's acclimated T crit (in green) under the same scenarios (CMIP6 ensemble median in black, individual ensemble members in red).T leaf-xx exceeds median T crit frequently after ~2050 (dotted vertical black line) under both SSPs suggesting conditions that maximise the potential for sustained damage to Larix sibirica's photosynthetic apparatus.

Fig. 4
Fig. 4 Projections of rising intensity and frequency of leaf temperature (T leaf ) in excess of the Larix sibirica's T crit .The percentage of ensemble members and the median number of days per year when T leaf exceeds different T crit thresholds under SSP5-8.5 (a, c) and SSP3-7.0(b, d).Thresholds shown are the minimum, 10th percentile, 25th percentile and median T crit values, assuming that T crit also acclimates to warming temperature.Results are computed across 56 ensemble members from 23 CMIP6 earth system models (ESMs).Despite the exceedance of T crit by T leaf generally being a lowprobability extreme event due to natural year-to-year climate variability, increases are projected both in the likelihood (percent of ensemble members) and magnitude (number of days/year) of T crit exceedances each year starting ~2050 across multiple models and model ensemble members.The solid line in a and b is a 3-yr running average.