Loss of vegetation functions during the Paleocene–Eocene Thermal Maximum

Domain: earth_science
Published: 2025-11-27
Editorial outcome: minor_revision

Abstract
The Paleocene–Eocene Thermal Maximum (PETM) around 56 million years ago was a 5–6 °C global warming event that lasted for approximately 200 kyr. A warming-induced loss and a 70–100 kyr lagged recovery of biospheric carbon stocks was suggested to have contributed to the long duration of the climate perturbation. Here, we use a trait-based, eco-evolutionary vegetation model to test whether the PETM warming exceeded the adaptation capacity of vegetation systems, impacting the efficiency of terrestrial organic carbon sequestration and silicate weathering. Combined model simulations and vegetation reconstructions using PETM palynofloras suggest that warming-induced migration and evolutionary adaptation of vegetation were insufficient to prevent a widespread loss of productivity. We conclude that global warming of the magnitude as during the PETM could exceed the response capacity of vegetation systems and cause a long-lasting decline in the efficiency of vegetation-mediated climate regulation mechanisms.

Introduction
During the onset of the Paleocene–Eocene Thermal Maximum (PETM), a geologically abrupt release of several thousand petagrams of isotopically light carbon, likely of mixed methanogenic, organic and volcanic origin1, increased atmospheric CO2 by more than 1000 ppm and global average temperatures by 5–6 °C1,2. With the estimated mass of carbon release and the magnitude of warming being comparable to future anthropogenic emission scenarios1, the PETM represents one of the most important geological analogues to present-day climate change and is key for understanding carbon cycling and biological dynamics in response to a period of severe global warming1,3.

Paleobotanical records from several locations around the globe suggest that the PETM induced a major disruption and reorganisation of terrestrial vegetation systems through plant range shifts and changes in species composition4–8, potentially impacting their climate-regulating functions9. Geochemical proxies of the atmosphere-ocean δ13C signature during the period (Fig. 1A, B) indicate that changes in vegetation-mediated organic carbon cycling and storage may have strongly shaped the PETM carbon dynamics. Similar to other large carbon reservoirs, such as methane hydrates10 or permafrost carbon11, organic carbon stored in vegetation and soils may have acted as a carbon cycle capacitor9,12: contributing to a global negative δ13C isotope excursion with the release of isotopically light carbon during the onset of the hyperthermal, remaining in a state of reduced carbon storage during an approximately 70–100 kyr period of low δ13C (the PETM body), and contributing to the carbon cycle and climate recovery as vegetation and soil carbon stocks regrew during the final stages of the event (the PETM recovery). High rates of δ13C recovery during the termination of the PETM carbon isotope excursions support that a recovery of the photosynthetic sequestration of isotopically light carbon may have contributed to the termination of the PETM9,12,13. However, the 70–100 kyr lag in the recovery of climate and isotopic composition suggests that the PETM warming shifted terrestrial ecosystems into a long-lasting state of reduced carbon storage, and that the time needed for vegetation systems to recover from the perturbation and to re-establish their carbon sequestration potential may have strongly influenced the duration of the hyperthermal.

Vegetation functioning in response to global warming depends on the ability of plants to avoid stress conditions through migration or to prevail in the new environment through adaptation and evolution14,15. For selected fossil species, large dispersal distances have been observed in response to the PETM warming5,16. However, it is unclear whether a widespread migratory response and the resulting reorganisation of vegetation structures6 were sufficient to attenuate the effects of the warming and to compensate for a potential loss of vegetation adaptation to the new environmental conditions. The evolutionary capacity of plants to adjust to new climates is expected to be limited in speed and magnitude, with climatic niches occupied by species remaining conserved on geologic timescales17–21. Slow evolution and a warming-induced adaptation lag could have reduced organic carbon sequestration during the PETM, and may have impacted additional vegetation-mediated carbon cycle feedback mechanisms, such as the biotic enhancement of CO2-consuming terrestrial silicate weathering reactions22–26.

Here, we assess how vegetation dispersal and evolution in response to the PETM affected the global carbon cycle and vegetation-mediated climate regulation. For this purpose, we developed a trait- and optimality-based biogeodynamical vegetation model, designed for continuous simulations of vegetation structure and functioning over geologic timescales. In addition to global simulations, for three sites in the mid- to high latitudes, we compare simulated vegetation trends with functional trait reconstructions using nearest living relatives of fossil species observed in palynofloral records, as well as available macrofossil observations. Using model simulations and vegetation reconstructions, we investigate the possibility of a disrupted terrestrial carbon cycle due to a warming-induced vegetation adaptation lag.

Results
Vegetation adaptation dynamics pace terrestrial carbon cycling
The modelled response of vegetation to the PETM warming affects the global dynamics of terrestrial carbon storage and silicate weathering (Fig. 1). Assuming a limited evolutionary capacity of vegetation to adapt functional and climatic traits to the environmental changes results in a reduction of terrestrial carbon storage following the PETM onset, potentially lasting as long as the body of the carbon cycle perturbation (70–100 kyr). A loss and lagged regrowth of the biospheric organic carbon storage is in agreement with geochemical records for the period, suggesting a release and reduced storage of biospheric carbon during the onset and body of the PETM, followed by an enhanced sequestration of isotopically light carbon during the termination of the carbon cycle and climatic perturbation9,12,13. A lagged organic carbon regrowth in the range of 1000 to 2000 Pg is modelled under scenarios of limited evolution. In contrast, when assuming fast trait and climatic niche adaptation or no limitation in climatic tolerances (Supplementary Fig. 8), a less severe initial drop in carbon sequestration is soon compensated by a rapid regrowth and overshooting in global biomass carbon compared with pre-PETM levels. The overshoot in biomass indicates that under a reduced or absent adaptation lag, the temperature, precipitation and atmospheric CO2 conditions during the PETM body would support globally increased rates of primary productivity. Such a vegetation–climate equilibrium, however, contrasts with fossil plant and soil records, which point towards a complete turnover of vegetation systems and diminished carbon storage following the PETM warming4–6,38.

The long duration of the recovery of carbon stocks is also observed in simulations considering the highest dispersal capacities inferred from the fossil record. Therefore, a migratory response may not have been sufficient to compensate for a loss of functional and climatic trait adaptation following the PETM environmental changes. Under a high dispersal scenario, carbon stocks are modelled to steadily recover from the perturbation, due to a continuous spatial redistribution of vegetation and a reduction of the adaptive lag. For a low dispersal scenario, there is a limited recovery of carbon stocks during the PETM body, followed by an acceleration during the PETM recovery period. The accelerated recovery under low dispersal occurs because modelled plants tend to remain close to the original habitat and rapidly exit the PETM stress state with the return of climatic conditions similar to those of the late Paleocene. Such carbon dynamics would be expected if plants managed to persist in local refugia throughout the climatic perturbation.

A vegetation adaptation lag additionally impacts carbon sequestration through a reduced vegetation-mediated enhancement of terrestrial silicate weathering. Under a scenario of slow vegetation recovery, the effects of a warmer and wetter climate on silicate weathering are temporarily offset by a reduced vegetation-mediated weathering enhancement following the PETM climatic perturbation. Due to the vegetation disruption, an overshoot in silicate weathering, and thus effective long-term climate regulation, may only occur with a lag of several tens of thousand of years.

Palynofloral vegetation reconstructions
We use palynofloral records of species composition and abundance spanning the PETM to reconstruct vegetation trait changes at the two mid-latitude sites Bighorn Basin6 and 22/11-N7, as well as the high latitude site ACEX (Arctic Coring Expedition)8 (Fig. 2A). Thereby, we assign the recorded fossil species an average trait value of plant height, leaf mass per area (LMA) and phenology (deciduous or evergreen), as observed in modern relatives of the same genus or family. Due to continuous evolution and adaptation, absolute trait values may have changed between the fossil species and their living relatives today. However, by assuming that the relative ranking in trait values remained conserved across clades (Supplementary Fig. 9), we aim to approximate major trends in vegetation structural changes across the PETM.

For both mid-latitude sites Bighorn Basin and 22/11-N, the palynofloral reconstruction indicate a shift towards smaller plants (Fig. 2), likely associated with a reduced biospheric carbon storage (Supplementary Fig. 10). For the Bighorn Basin, the reconstructed reduction in vegetation height coincides with a reduced abundance of deciduous trees and an increase in average LMA (smaller and thicker leaves). The trait shifts are caused by a dominance of tropical to subtropical evergreen palms (i.e., Serenoa, Pseudophoenix) and ferns (i.e., Pityrogramma, Acrostichum, Cibotium, Salvinia, Anaemia) during the PETM body, while tall-growing trees of a reduced average LMA and likely more frequent deciduousness (i.e., Metasequoia, Platanaceae, Betulaceae, Juglandaceae) only reappear during the PETM recovery period. The vegetation shifts in the Bighorn Basin suggest a dominance of more conservative, but likely more robust plant strategies concerning potential heat and drought conditions during the PETM body39. Similarly, at the site 22/11-N, evergreen palms and ferns (i.e., Schizaeaceae, Polypodiaceae) and small wetland floras (i.e., Typhaceae) dominate the playnofloras of the PETM onset and body, resulting in a reduction of the average vegetation height. Mixed conifer-broadleafed vegetation (i.e., Pinaceae, Metasequoia, Betulaceae, Juglandaceae) reappears during the PETM recovery period. Compared to the Bighorn Basin, a less drastic height shift is reconstructed due to the persistence of some conifer and broadleaved trees throughout the perturbation.

The palynofloral reconstructions at both mid-latitude sites are best reproduced in the model simulations that consider a limited evolutionary capacity of vegetation to adapt functional and climatic traits following the PETM onset, resulting in a primarily dispersal-driven vegetation response. The loss of local adaptation results in reduced productivity, carbon sequestration and height growth during the body of the PETM. The combination of simulations and palynofloral records therefore support the possibility of a warming-induced adaptation lag, resulting in the dominance of robust, but more conservative plant strategies following the PETM onset, until more acquisitive strategies (i.e., tall, deciduous trees of low LMA) slowly return with time. Importantly, in contrast to the palynofloral reconstructions, simulations that do not consider an adaptation lag indicate a theoretically possible increase in vegetation height and carbon sequestration potential at both the Bighorn Basin and 22/11-N sites after the onset of the PETM warming. A loss of vegetation adaptation may therefore have offset environmental conditions in favour of increased primary productivity, causing more robust but less productive plant strategies to dominate the local vegetation. An additional indication of a limited evolutionary trait and climate adaptation capacity is the high degree of species turnover observed at the sites, with only a few species persisting throughout the climate excursion by either enduring the new climate or adjusting to it.

A mismatch is observed between modelled and reconstructed LMA for the 22/11-N site following the PETM onset. The model predicts an initial shift towards increased LMA that is not observed in the data, whereas both reconstructions and the slow adaptation model predict a subsequent decreasing trend during the PETM body. The mismatch in the initial shift is likely caused by the model not explicitly representing the herbaceous understory vegetation that in terms of abundance, dominates the PETM-body palynoflora of 22/11-N. Instead, the optimality-driven model only captures the plant growth form that is most competitive regarding carbon sequestration and height growth. The model predicts the presence of 10–15 m evergreen trees of high LMA during the PETM body. This is in the range of growth forms reconstructed for the 22/11-N site from the palynofloral records, although at a reduced abundance (Supplementary Fig. 12). In both cases, the reconstructed and modelled reduction of height growth during the PETM warming suggests a reduced carbon sequestration at the site, supporting an adaptive lag and disrupted vegetation state.

We observe contrasting patterns for the high-latitude site ACEX, where simulations and palynofloral reconstructions suggest an increase or fast recovery of vegetation height and biomass carbon storage following the onset and during the body of the PETM. Subsequently, a decrease in the vegetation height is observed during the termination of the hyperthermal and the associated climate cooling during the PETM recovery period. The trait trends are driven by a replacement of pre-PETM mixed conifer-broadleaved forests with abundant needle-leaved taxa (i.e., Pinus, Picea) by broad-leaved swamp forests (i.e., Juglandaceae), including some rare subtropical evergreen palms (i.e., Arecaceae) and deciduous tropical trees (i.e., Bombacoideae) during the body of the climate perturbation. These reconstructed and modelled trends indicate that increased temperatures (Fig. 2B), longer growing seasons, sufficient precipitation, and increased atmospheric CO2 levels may have enabled a fast recovery or prevented negative warming effects for vegetation systems in the high latitudes.

Comparison with macrofossil observations
The palynofloral trait reconstructions rely on trait conservatism through time and only capture trait changes driven by turnover in species composition. A direct comparison of model simulations with paleo-functional trait observations, and an assessment of potential evolution-driven changes in the trait ranking among the occurring species (Supplementary Fig. 9), is limited by the availability of preserved macrofossils. The Bighorn Basin is the only locality with a stratigraphically dense record of plant macrofossils throughout the PETM40. Macrofossils from the Bighorn Basin generally support the vegetation dynamics derived from the palynofloral record, indicating a replacement of mixed conifer-deciduous temperate forest taxa (e.g., Platanaceae, Betulaceae, Cupressaceae) with tropical and subtropical taxa (e.g., Fabaceae, Arecaceae, Hernandiaceae), likely resulting in a more open vegetation and an increased abundance of heat- and aridity-adapted species4,39,40. No direct estimates of height and biomass are available, but organic carbon measurements from terrestrial sediments and fossil soils indicate a major reduction in organic carbon contents during the PETM41,42 (Supplementary Fig. 15). Depleted soils and reduced carbon inputs from plant productivity are consistent with the reduced productivity and height growth derived in the palynofloral reconstruction and the slow evolution model (Fig. 2). Petiole width-based estimates of LMA for dicot angiosperms in the Bighorn Basin between the late Paleocene and early Eocene show no significant within-species LMA variations43,44, which supports limited evolution-driven trait changes on PETM timescales. Further, similar as in the palynofloral reconstruction, the PETM indicates a marked increase in the abundance of high LMA species (i.e., Fabaceae) (Supplementary Fig. 15).

For the high latitudes, the temporal dynamics reconstructed and modelled for the site ACEX during the PETM resemble the latitudinal vegetation gradient observed in macrofossils from the late Paleocene and early Eocene Canadian Arctic. Despite the extreme photoperiodicity in these latitudes, macrofossils indicate forests composed of conifer and broadleaved (low LMA) deciduous trees, as well as some tropical palms in the lower Arctic latitudes45,46. Leaf physiognomic analyses indicate growth conditions similar to today’s temperate forests or temperate rainforests45. Fossil tree stumps from the late Paleocene show diameters that indicate tree heights well beyond 20 m, in line with the height growth potential derived from the palynofloral reconstruction47. While there is no macrofossil assemblage available that confidently capture the onset, body and recovery phases of the PETM, there is currently no Canadian Arctic late Paleocene to early Eocene macrofossil site indicating a pronounced loss of productivity or diversity during the period45,48. This is supported by limited changes observed in the organic geochemical composition of coals deposited in the Canadian Arctic across the PETM49, and agrees with the stable vegetation reconstruction derived from the palynofloral species composition at ACEX.

Warming scenarios and CO2 fertilisation
Modelled vegetation functioning during the PETM depends on the magnitude of experienced surface warming. We tested the vegetation response to different land surface warming scenarios, from an average of 1 °C to 8 °C (PETM body CO2 concentrations between 800 and 1590 ppm), with the warmer endmember corresponding to estimates for the PETM2 (Fig. 3A, B). From our simulations, we estimate a threshold around 4 °C, beyond which a widespread migratory response vegetation is not able to compensate negative warming impacts, resulting in an adaptation lag and a reduction of vegetation productivity levels. Warming beyond this threshold results in non-analogous climatic environments, making vegetation-mediated carbon sequestration primarily dependent on the capacity of plants to adapt traits through evolution and to occupy climatic niches well outside their usual habitat. Below this warming threshold, no adaptation lag is observed, and increased atmospheric CO2 concentrations and warmer climates promote increased levels of primary productivity.

We further find that the temporal dynamics of terrestrial carbon storage depend on the vegetation’s ability to use increased atmospheric CO2 levels for photosynthesis during the PETM climatic perturbation (Fig. 3C–E). The CO2 fertilisation effect reduces the water loss per unit of carbon taken up for photosynthesis, increasing the vegetation’s resilience to drought and heat. Our simulations indicate that any limitation in the CO2 fertilisation effect strongly reduces the vegetation’s carbon assimilation potential, and prolongs the lag in the recovery of terrestrial carbon stocks and productivity. A limitation in the CO2 fertilisation effect could have been an additional consequence of a loss of local adaptation following the PETM climatic perturbation, resulting in the inability of plants to efficiently acquire and utilise local water and nutrient resources for photosynthesis. The importance of CO2 fertilisation to maintain productivity and biomass is particularly pronounced in the mid-latitudes.

Methods
Vegetation model
Vegetation modelling was conducted using the newly developed TREED model (TRait Ecology and Evolution model over Deep time, version 0.1). In addition to the summary of the model provided here, an extended description of all model equations and parameters is given in the Supplementary Methods (SM). TREED is a trait- and optimality-based vegetation model designed to conduct continuous vegetation simulations over geologic timescales. Based on concepts of eco-evolutionary optimality29,30, the model considers variable, continuous and evolving functional and climatic trait combinations of modelled plants over time. Following the scheme in Supplementary Fig. 16, the model consists of four modules executed in succession.

Initialisation, allometric and functional relationships. The model is initialised with one reference plant per terrestrial grid cell, characterised by six key traits that describe the vegetation structure and functioning of the local vegetation. The key traits include: plant height (H), mass of the leaf carbon pool (Cleaf), phenology (deciduous or evergreen), fineroot to leaf carbon ratio (r:s), leaf longevity (all), and climatic niche. The climatic niche describes the climatic environment to which the modelled plant is best adapted, characterised by the temperature of the coldest month (Tmin,veg), the temperature of the warmest month (Tmax,veg), and the mean annual temperature (Tmean,veg) of this environment. Large differences between local climatic conditions and a modelled plant’s climatic niche, for example, after climatic changes or migration, can result in a reduction in productivity, as described below. At model initialisation, the climatic niche of the modelled reference plant in a grid cell is derived from local climatic conditions, assuming an optimal adaptation. Using the mentioned key traits and allometric and functional relationships, several other plant characteristics relevant to vegetation structure and functioning are derived. One key trait relationship is the dependency between a modelled plant’s specific leaf area (SLA) and leaf longevity (all), according to (all model parameters and variables, including units, are listed in Supplementary Tables 1 and 2):

SLA = (2 * 10^-4) * 1 * 10^(2.25 - 0.5 * log(all * 12)) / DMc (1)

This functional relationship represents a generalised form of the relationship used for needle-leaved and broadleaved plants in the LPJml4 dynamic vegetation model28. It represents a carbon-economic trade-off between a high SLA, promoting a large leaf area and thus a high carbon capture potential, and high annual leaf building costs due to a low leaf longevity33. DMc is an assumed averaged dry matter carbon content of 0.47 g C per g of dry matter, following28. SLA and carbon assimilation potential are related through the leaf area index (LAI), which ultimately determines the potential to capture light for photosynthesis:

LAI = Cleaf * SLA (2) / CA

where CA is the “crown area", representing the surface area occupied by a modelled plant and that scales with plant height (SM).

Several additional allometric relationships are used to describe plant geometry and the carbon partitioning between different plant carbon pools (Csapwood, Cheartwood, Cfineroot, Ccoarseroot, Cleaf) (SM). Given that there are no pre-defined plant functional types, the allometric relationships are generalised for all modelled plants. All parameters related to geometry and carbon partitioning are calibrated using present-day vegetation height and biomass data.

Carbon balance. Based on a reference plant’s traits, the carbon balance module estimates photosynthetic carbon assimilation, maintenance and growth respiration, as well as average carbon investments for tissue turnover at a location. Photosynthesis is modelled based on the Farquhar photosynthesis model, considering generalisations for global modelling applications27,28,73–76. Based on ref. 28, the absorbed photosynthetically active radiation (APAR) in each terrestrial grid cell depends on the incoming shortwave radiation and the light capture potential of the local reference plant, according to:

APAR = PAR * (1 - e^(-k * LAI)) * αleaf:stand * Φnichestress (3)

where half of the downwelling surface shortwave radiation (RSDS) is assumed to be photosynthetically active radiation (PAR = 0.5 * RSDS) and αleaf:stand is a leaf to stand level photosynthesis scaling parameter. k is a light extinction coefficient, which is 0.6 for deciduous plants with all < 1, 0.4 for plants with all > 4, and 0.5 in all other cases. Φnichestress is an introduced stress factor that represents a potentially reduced carbon capture potential for plants in climatic environments to which they are not adapted. It is calculated as:

Φnichestress = min(Φheat, Φcold, Φmean) (4)

with

Φheat = exp(-knichebreadth * (Tmax,local - Tmax,veg)^2) (5)

Φcold = exp(-knichebreadth * (Tmin,local - Tmin,veg)^2) (6)

Φmean = exp(-knichebreadth * (Tmean,local - Tmean,veg)^2) (7)

where k describes the impact of a deviation of the local warmest month temperature (Tmax,local), local coldest month temperature (Tmin,local), and local mean annual temperature (Tmean,local) from the plant’s climatic niche (Tmax,veg, Tmin,veg and Tmean,veg, respectively). Three k values are tested, i.e., 0, 0.03 and 0.10, corresponding to no, intermediate and strong limitation of primary productivity in new climatic environments, respectively (Supplementary Fig. 8). The default intermediate niche breadth is derived from the temperature variations observed within present-day biomes (i.e., regions with similar vegetation types) (Supplementary Fig. 1) and considers a maximum 10 °C deviation in local mean, warmest month or coldest month temperature from the plant’s climate niche before a loss of productivity occurs. Based on APAR, photosynthesis is calculated as the minimum of light-limited and Rubisco-limited photosynthesis (SM). The model has the form of a “light-use efficiency” model, with PAR, temperature, day length, water availability and canopy conductance as the major modulating factors.

The rate of carbon assimilation is subject to a closed water balance, where the water loss from vegetation during photosynthesis (actual evapotranspiration, AET) cannot exceed the water supplied by precipitation. Following27, AET is calculated as:

AET = min(Esupply, Edemand) (8)

where Esupply is the monthly average precipitation rate, and a monthly average Edemand is calculated following77:

Edemand = Eeq * (1 + αm * gc / gmc) (9)

where αm (= 1.391) is a Priestley-Taylor coefficient, gm (= 3.26 mm/s) is a scaling factor28, gc is the canopy conductance, and Eeq the equilibrium evapotranspiration rate (SM). Canopy conductance and photosynthesis are related by27,28:

gc = 1.6 * A / (p * (1 - λ)) + gdt (10)

where A is the daily net daytime photosynthesis (SM), gdt is the minimum canopy conductance that occurs due to non-photosynthesis-related processes, and pa is the ambient partial pressure of CO2. λ describes the leaf–atmosphere water and carbon exchange and equals λmax (0.8) under non-water-limited conditions. If Edemand > Esupply, λ is adjusted to λ < 0.8. In cases of water limitation, λ is solved for to fulfil Eq. (10) and to obtain a closed water balance (Edemand = Esupply). The adjustment of λ results in an adjustment of the photosynthetic rate and the canopy conductance accordingly.

Maintenance respiration rates for metabolically active carbon pools (Cleaf, Csapwood, Cfineroot) depend on temperature and tissue-specific C:N ratios for sapwood and fineroot carbon, and on the maximum Rubisco capacity for leaves (SM). A fixed fraction of the modelled plant’s net carbon assimilation is assumed to go into growth respiration27,28, resulting in:

NPP = (1 - rgr) * (GPP - Rleaf - Rsapwood - Rfineroot) (11)

where rgr is the fraction that goes into growth respiration and GPP (gross primary productivity) is the annual sum of the monthly average rates of leaf, sapwood and fineroot respiration.

The model does not explicitly resolve the establishment, growth and mortality of individual plants. To evaluate whether a certain vegetation form is sustainable in a location, an annual average carbon balance is calculated, considering the average carbon costs of tissue turnover derived from the modelled reference plant’s characteristics and carbon pools. The carbon balance ensures that the carbon capture potential of the reference plant is large enough to build and maintain the modelled plant characteristics. For sapwood, heartwood and coarseroot carbon pools, the carbon costs of tissue turnover (τsapwood, τheartwood, τcoarseroot, respectively) are calculated based on tissue-specific turnover times (f) and heat damage (HD) and frost damage (FD) indices, accounting for higher carbon costs of turnover in extreme environments (maximum doubling; SM), resulting in:

τsapwood = Csapwood * fsapwood * (1 + HD + FD) (12)

τheartwood = Cheartwood * fheartwood * (1 + HD + FD) (13)

τcoarseroot = Ccoarseroot * fcoarseroot * (1 + HD + FD) (14)

The carbon costs of leaf turnover depend on leaf longevity:

τleaf = Cleaf * (1 / all) (15)

Further, for leaf longevities longer than a year, the carbon costs of fineroot and leaf turnover are assumed to be proportional28:

τfineroot = Cfineroot * (1 / all) (16)

Subsequently, the annual carbon balance can be calculated as:

NCG = NPP - (τleaf + τsapwood + τheartwood + τcoarseroot + τfineroot) * 1 / CA (17)

where NCG is the net carbon gain, i.e., the carbon available from photosynthesis after accounting for the carbon costs of maintenance and growth, respiration and tissue turnover. This carbon is available to the plant for reproduction, defence, nutrient acquisition and other processes relevant to survival. NCG is thus considered a fitness measure to compare the suitability of plant characteristics in a given environment.

Trait optimisation. Over time, environmental changes result in selection pressures that drive changes in vegetation structure and functioning. Following optimality principles29,30, it is assumed that natural selection results in trait combinations evolving towards increasingly optimal adaptation to the local environment. To imitate this selection process, an optimisation algorithm is employed at every timestep and for every terrestrial grid cell to predict the trait combination that maximises NCG (considering radiation, temperature and precipitation as environmental inputs). For this, a differential evolution algorithm for multi-parameter optimisation problems (DEoptim) is employed78. The optimisation predicts the optimum Cleaf, phenology and all. In environments with more than 3 months with average temperatures below 3 °C (considered the necessary condition for a deciduous phenology), the optimisation procedure evaluates whether deciduousness results in a larger NCG compared with an evergreen phenology. In this model step, the maximum potential plant height for a given trait combination is also evaluated, defined as the height at which the ratio of NCG to NPP is approximately 0.2, following the assumption that this is the minimum NCG necessary to ensure reproductive success.

Adaptation and dispersal. At every timestep of the model (10 kyr), the traits of the modelled plants evolve toward the predicted local optimum trait combination from the optimisation procedure. The rate of trait adaptation per timestep is a user-defined adaptation rate (α), which is varied between model simulations in the range from 0.01 to 1 to estimate the effect of a low to high ability of vegetation to adapt to new environments through trait evolution. Thereby, an adaptation rate of 1 indicates a trait change per timestep that amounts to 100% of the difference between the predicted optimum and the current trait value, whereas an α of 0.01 represents an evolutionary trait change amounting to 1% of the difference. In the present model, the traits Cleaf, all and climatic niche are subject to α. The climatic niche is assumed to evolve towards local conditions, resulting in a continuously more specialised adaptation to the local environment the longer a reference plant persists in a locality:

Cleaf,new = Cleaf,old + N(α, 0.05) * (Cleaf,target - Cleaf,old) (18)

all,new = all,old + N(α, 0.05) * (all,target - all,old) (19)

Tmin,plant,new = Tmin,plant,old + N(α, 0.05) * (Tcoldestmonth,local - Tmin,plant,old) (20)

Tmax,plant,new = Tmax,plant,old + N(α, 0.05) * (Twarmestmonth,local - Tmax,plant,old) (21)

Tmean,plant,new = Tmean,plant,old + N(α, 0.05) * (Tmean,local - Tmean,plant,old) (22)

Where “target” denotes the trait value predicted to be optimal from the optimisation procedure, and “local” is the local climatic variable. To allow a limited degree stochasticity in the evolutionary adaptation in space and time N(α, 0.05) indicates that the rate of evolution is drawn from a normal distribution with mean α as defined for the simulation run and a standard deviation of 0.05. The height trait is considered to be fully dynamic and set to the evaluated maximum potential height with a plant’s current trait combination, accounting for the dynamic response of plant heights to environmental conditions. However, given that Cleaf is subject to α, there is some degree of conservatism in the growth form between time steps. In locations where a deciduous phenology is assessed to result in a favourable NCG, and where all,new is ≤ 1 year, the modelled plant is considered deciduous.

In addition to trait adaptation, vegetation can respond to environmental changes through dispersal. At every timestep, the modelled plants can disperse into the surrounding grid cells within a user-defined range. By testing variable dispersal capacities of the vegetation, the importance of trait adaptation vs dispersal and redistribution dynamics is tested. Dispersal is modelled as a stochastic process, accounting for the large range of factors involved in plant dispersal but not explicitly represented here (e.g., dispersal vectors, weather, soil properties)79. At every timestep and for every modelled plant, a new dispersal distance is drawn from a dispersal kernel, defined by a Weibull distribution with a shape parameter value of 1.75 and a user-defined scale parameter. In this study, the dispersal scale is varied between simulations in the range from 200 to 1100 km per 10 kyr, covering a large range in possible dispersal capacities as inferred from plant fossils for the Paleocene–Eocene Thermal Maximum (PETM)5. The shape of the dispersal kernel accounts for the general right-skewed distribution of dispersal distances, with a high frequency of low dispersal distances and a low frequency of long-distance dispersals67. During dispersal, offspring with trait combinations identical to those of the reference plant migrate into the grid cells within the distance of the drawn dispersal value.

If several reference plants with a positive NCG occur in a location after dispersal, a competition function selects the surviving reference plant with the best trait combination. It is assumed that the reference plant with the highest NCG and plant height will dominate the local vegetation, considering carbon gain and success in light competition as the main fitness measures. To account for the combined height and NCG competition during the selection, NCG is scaled with a height penalty function that depends on the average height of the competing reference plants34.

To model trait evolution, dispersal and competition dynamics, the described model functions are implemented in the “general engine for eco-evolutionary simulations” framework (gen3sis)31. The framework enables efficient handling of the described eco-evolutionary steps while considering a spatially-explicit dynamic landscape with temporally variable climatic input data.

The choice of the timestepping in the TREED model is variable and user-defined. In this study, a timestep of 10 kyr was employed to resolve the abruptness of the PETM onset climate warming, while still being able to cover the full range of the PETM of 200 kyr, and to test a large range of dispersal, evolution, warming and CO2 fertilisation scenarios. Thus, in the presented implementation, at every 10 kyr time step, the model takes multiyear monthly average climatologies as inputs and calculates average vegetation structures and carbon fluxes, but does not resolve variability at any other temporal resolution in between 10 kyr time steps and apart from the seasonal variability derived from the monthly climate inputs.

Soil organic carbon. Soil carbon stocks are based on a reduced-complexity soil carbon model derived from35. In this model, a one-dimensional soil layer at steady state is considered, where inputs from NPP are equal to heterotrophic respiration. It is assumed that soil heterotrophic respiration is directly proportional to the soil carbon pool and modulated by temperature according to a Q10 function with a baseline temperature of 15 °C80:

Csoil = NPP / (k * Q10^((T-15)/10)) (23)

with a uniform k of 1 year^-1 and a Q10 of 1.75.

Climate data
The TREED vegetation model is forced with monthly output from climate simulations using the Community Earth System Model version 1.2 (CESM1.2)6. For the CESM1.2 climate simulations, a high-resolution fixed sea surface temperature configuration with coupling between the atmosphere component, CAM5.3 (Community Atmosphere Model version 5.3), the land component CLM4 (Community Land Model version 4), and the River Transport Model (RTM)81–83 was applied. Sea surface temperatures were derived from lower-resolution (≈2°) fully-coupled companion simulations, run to equilibrium (≈1800 years). CESM1.2 model boundary conditions were taken from Deep-MIP forcing datasets84,85. Two climate simulations are used as the basis for the vegetation modelling, one with a pre-PETM atmospheric CO2 concentration of 680 ppm and one with a PETM CO2 of 1590 ppm. The simulations were carried out assuming a neutral orbital configuration (0 eccentricity, obliquity of 23.5°).

For the continuous TREED vegetation model simulations, climate conditions before the PETM onset are derived from the CESM1.2 pre-PETM simulation at 680 ppm atmospheric CO2 and are kept constant during a 50 kyr period. After this period, the TREED model is forced with the climate conditions of the PETM body (1590 ppm CO2), with a linear transition of climate outputs (duration of 10 kyr). Then, assuming a PETM body duration of 70 kyr, CO2 concentration and climate are kept constant at 1590 ppm CO2. Finally, during the PETM recovery period, the CO2 concentration is assumed to linearly decrease to 780 ppm CO2 until 150 kyr after the PETM onset, and to pre-PETM conditions of 680 ppm CO2 until 200 kyr after the PETM onset. For the transition climates, the approach from other deep-time biogeochemical models86,87 is followed, and climatic fields are interpolated from the available climate model runs. The same procedure is applied for testing the effects of different land surface warming scenarios on vegetation dynamics (i.e., Fig. 3). The PETM CO2 forcing and the resulting climatic fields are illustrated in Supplementary Fig. 7.

Palynofloral trait reconstructions
The palynofloral records compiled in ref. 6 are used for the PETM vegetation trait reconstructions. Only records with a degree of age control that allow differentiation of “PETM body” and “PETM recovery period” samples are considered, in agreement with the carbon isotope records shown in Fig. 1. To derive vegetation traits from the fossil spore and pollen record, the nearest living relative (NLR) approach is used. Instead of assuming trait conservatism in the climatic ranges between fossil species and NLRs, as in the traditional NLR approach6, we assume conservatism in plant traits. We consider the two approaches similar in their underlying assumption, given that plant traits ultimately determine the ability of species to persist in a given environment. In the trait estimation procedure, each fossil species is assigned an average trait value derived from observed trait values of their NLR under present-day conditions. The NLRs for the fossil species observed in the spore and pollen data are derived from ref. 6. Present-day trait values are derived from the communal plant trait database TRY88, with a full reference list of the considered trait measurements given in the Supplementary References. The traits plant height, leaf mass per area (=1/SLA), and phenology are considered. Regarding phenology observations in the TRY database, a value of 0 is assigned to observations of deciduousness, and 1 assigned to observations of an evergreen habit. From the available trait data, average values are calculated on the family and genus level and assigned to the fossil species of the same family or genus. With every observation in the TRY database being given the same weight, the average will approximate the most common trait value within a genus or family. While not representing the full range of variability observed in natural systems, or within a family or genus, the generalised approach aims to approximate the major ranking of trait values among the involved taxa (Supplementary Fig. 9). If a fossil species is assigned several possible nearest living relatives, an average value across all observations from the plausible genera and families is assigned. Genus-based trait assignments may give more specific trait values than family ones, due to a higher degree of trait conservatism at a lower taxonomic level. However, not all fossil species have living relatives at the genus level, resulting in a low coverage of trait values for the observed fossil species. This is why both genus-based and family-based reconstructions are presented. To give a temporal evolution of the average trait combination at the different localities, the assigned trait values are weighted by the abundance of the fossil species in the sample. An illustration of the full range of trait values in the different samples and their relative abundances is given in Supplementary Figs. 11, 12 and 13.

Silicate weathering estimation
Silicate weathering rates throughout the PETM are estimated considering the CESM climate inputs, the topographic reconstruction of ref. 84 and vegetation biomass estimates from the TREED vegetation model. We apply a weathering equation commonly used in deep-time biogeochemical models86,89:

ωsilw = χm * ϵ * (1 - exp(-K * e^(Ea * (T0 - T) / (R * T * T0)))) * (1 - e^(-kw * q)) * (z / (σ + z)) * (1 + fBM) (24)

where ωsilw is the cation flux released during silicate weathering and ϵ is the local erosion, calculated as a function of topographic slope (s), runoff (q), temperature (T) and a scale parameter (k; calibrated to reproduce erosion rates similar to present-day levels) following90:

ϵ = ke * q^0.31 * s * max(T, 2) (25)

In Eq. (24), T is the local temperature, Xm the cation abundance in bedrock (0.1), z the regolith thickness (10 m), Ea the apparent activation energy of silicate weathering (20 KJ mol^-1), R the ideal gas constant, and T0 the standard temperature. K, kw and σ are calibration constants (6 * 10^-5, 10^-3 and -0.1, respectively, adopted from ref. 86). Finally, fBM is a vegetation-dependent weathering enhancement factor. Based on the “carbon-energy flux” hypothesis24, weathering enhancement is assumed to scale with carbon and nutrient fluxes, which intensify with biomass, resulting in a weathering enhancement formulation derived from ref. 36:

fBM = [1 - min(BMnorm, 1)] * NOPLANT + BMnorm (26)

where BMnorm is a normalised biomass indicator between 0 and 1. The biomass indicator linearly scales with the biomass density derived from the TREED model, with 0 indicating no biomass, and 1 indicating biomass densities as observed in today’s tropical forests (140 Mg C per ha, or around 300 Mg dry biomass per ha)91. NOPLANT describes the background weathering intensity in the absence of vegetation. In agreement with laboratory and field studies, a maximum six-fold weathering enhancement by vegetation is assumed (NOPLANT=1/6, 22, 23, 25). The CO2 consumption by silicate weathering (Fsilw) is assumed to be proportional to the cation release ωsilw 86, with 1 * 10^13 mol C year^-1 as the baseline carbon consumption (kref; in the range of present-day carbon consumption rates through silicate weathering) and ωref as the baseline cation release flux under the initial climatic conditions:

Fsilw = kref * (ωsilw / ωref) (27)

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Rogger et al. present results from a new eco-evolutionary vegetation model (TREED), driven by boundary conditions from climate model simulations using the Community Earth System Model of the pre-PETM and PETM, in order to constrain the vegetation response across this past global warming event. The authors evaluate the model against proxy reconstruction of vegetation response at three sites (two in the Northern mid-latitudes and one in the Arctic) based on pollen assemblages and using nearest living relatives to reconstruct modeled vegetation traits. The experimental framework evaluates the relative role of adaptation speed and dispersal as controls on the vegetation response and contrasts cases with and without CO2 fertilization and different degrees of warming. Based on the response of vegetation, the model further constrains the impact on terrestrial carbon storage (by combining the vegetation model with a simple representation of soil carbon storage) and the impact on silicate weathering (by adding a scaling of the weathering function based on biomass density). Combined, the model sensitivity tests and comparison with vegetation reconstructions from the three sites suggest that PETM warming exceeded a threshold for the ability of vegetation to adapt or migrate, thus leading to a temporary loss in productivity and hence total land carbon storage as well as a counteracting influence on the warming-induced increase in silicate weathering. The modeled vegetation response is consistent with both a long-duration of the so-called ‘PETM body’ and abrupt recovery.\n\nI think this is a great study and very useful modeling contribution, with relevance both for the PETM and for the long-term future response of vegetation to anthropogenic warming. The manuscript is well written, the model description is thorough and the presentation of results is easy to follow, both in the text and the figures. At the same time, there is a lot included in the manuscript, so I have a few suggestions, which may help contextualize the results. Overall, I think this paper is really well done and worthy of publication in Nature Communications.\n\nFirst, in a few places in the introduction (lines 46, 52, and 56) the text doesn’t clearly discriminate between biotic carbon stocks versus terrestrial carbon stocks specifically. Slight changes to the wording could clarify that existing d13C records do not uniquely implicate terrestrial organic carbon as a capacitor, as opposed to other organic carbon sources (marine organic carbon, hydrates, permafrost carbon, etc.). Changes in terrestrial organic carbon storage (vegetation and shallow soils, excluding permafrost) are one possibility, and this paper is evaluating the hypothesis that the terrestrial ecosystem response can account for the highlighted aspects of the global d13C records.\n\nNext, the Methods section answered most of the model-specific questions that I had while reading the text, but I think it would be helpful to include descriptions of the following in the main text (‘Trait based vegetation model’ section) in a third paragraph before Results: 1) a brief description of how the vegetation model interacts with soils - e.g. the sentence in the Figure 1 caption that soil carbon is a function of net primary productivity and temperature-dependent heterotrophic respiration, 2) a brief description of how silicate weathering is calculated - e.g. move forward the description on Line 160 and clarification of whether soils and weathering are part of the vegetation model or computed offline, 3) a brief description of the experimental framework, e.g. how the time-dependent changes were constrained (50 kyr using pre-PETM constraints, 70 kyr of PETM constraints, then a linear CO2 decrease over 130 kyr with linear interpolation of climatic fields), and 4) a brief summary of the sensitivity tests - varying the adaptation and dispersal capacity and including CO2 fertilization and different degrees of warming. Specifically, how are different amounts of warming imposed? What is the amount of warming that is consistent with the CESM simulations (pre-PETM to PETM)? Is scaling applied to both spatial patterns in temperature change as well as the seasonal temperature range? These additional descriptions of both the model and the experimental framework would allow much more clarity in interpreting the subsequently presented results.\n\nA few more specific model questions (but that could go into the Methods):\nLine 99: are modern orbits assumed in calculating seasonality in both the pre-PETM and PETM simulations with CESM?\nLine 523: what is the rationale behind selecting a 10 kyr model timestep?\nLine 572: what is the temporal variation of climate input data? Monthly?\n\nIn the discussion, it would be helpful to add some text to clarify the lack of feedback between the imposed climate boundary conditions and vegetation response and consider implications of that assumption. For instance, a modeled drop in biomass should generate not only a CO2 feedback via carbon release but also have a climate influence via changes in albedo, evapotranspiration, etc., and these cannot be represented with this experimental framework. Such feedback could potentially influence the timescale of the vegetation-mediated climate response.\n\nFurther, how have the results been shaped by assumptions in the modeling framework - specifically, the imposition of the PETM climate for 70 kyr? What if, instead, the experiment imposed a linear recovery in CO2 and temperature without a lag and more closely approximating the modeled response to a single pulse release of CO2 as has been simulated in models of the ocean-atmosphere (e.g. cGENIE, LOSCAR)? Would there still be a prolonged reduction in productivity, which in turn might explain why the actual PETM warmth persisted over a long duration? Given the number of sensitivity tests already included in the study - evaluating trait adaptation capacity, dispersal capacity, CO2 fertilization and degree of warming - I am not suggesting that it would necessarily add clarity to evaluate different assumptions about the climate boundary conditions, but it would be useful to discuss how the assumptions about these boundary conditions have impacted the results.\n\nFinally, in making comparisons between results for the PETM and the long-term future, do the authors have a sense whether the response is related to the relative or total amount of warming? Is 4 deg C a threshold because of an absolute temperature reached, or would 4 degrees of warming from a lower baseline have similar results, particularly because of the climatic niche constraint?", "reviewer_2": "Review of: “Loss of vegetation functions during the Paleocene-Eocene thermal maximum” by Rogger, J., Korasidis, V. A., Bowen, G. J., Shields, C. A., Gerya, T. V., and Pellissier, L., manuscript submitted to Nature Communications\nI read with much interest the manuscript by Rogger et al. It is a quite interesting study with significant novelty that deserves to be published in Nature Communications.\n\nThe authors present a new vegetation model with the originality that it is not based on pre-defined plant functional types, but on the optimization of plant traits to the local climate. This offers the advantage of plant traits varying in a continuous way in the modelled vegetation. Furthermore, the model can take into account adaptation lags under changing climate, as well as dispersal from nearby pixels. This is a rather unique tool to study the response of natural vegetation to palaeoclimate changes.\n\nThe model is tested over the period of Paleocene-Eocene Thermal Maximum (PETM). The PETM warming was likely too fast and too intense to allow the adaptation of the terrestrial vegetation. For global land surface warmings larger than 4°C, the model predicts a strong drop of net primary productivity and biomass of the terrestrial biosphere, that presumably impacted the carbon cycle and silicate weathering, implying a positive feedback on atmospheric CO2. This is an important result, which may have implications in the framework of the present climate warming.\n\nThe current version of the paper is already a revised version. I had not the opportunity to comment on the original version. Nevertheless, I checked that the authors responded completely to all reviewers’ comments. Their response looks fine to me.\nOne of the main points raised was linked to the reconstructions of palaeo plant traits from palynofloral records. Indeed, plant traits are reconstructed at the genus or family level (which is the best that can be done from palynoflora) and assumes that the average genus or family traits can be obtained from the modern TRY database, i.e., that these average traits at genus or family levels have not changed from the Paleocene-Eocene to the present. This may look inconsistent with the idea (and the model) that traits optimize in response to local climate. It is well-know that traits can vary strongly between individuals within a genus, or even within a single species. However, the current knowledge does not allow to make precise reconstruction of past trait values through times, especially for the PETM where macrofloral data may not provide the necessary resolution to assess the evolution through the PETM. Anyway, the authors now reconstruct LMA trait from macrofloral data at one of their sites (Bighorn site) where such a reconstruction is possible. The trends are consistent with the results obtained from the palynoflora. But, more importantly, the authors now state clearly their hypothesis, especially that the aim is not to reconstruct absolute values of traits (which is probably not possible), but relative changes over the PETM period. The hypothesis is then that the relative differences between genera or families have not changed since the Paleocene. This is not necessarily true, but it is fine provided that the hypothesis is clearly stated.\n\nI have noted a few other problems that I am summarizing below. The paper can be accepted after these small changes are taken into account.\n\nThe most significant problem is related to the validation of the model on the present (Supplementary figure 2). This is important, since the model is new and rests on a slightly different methodology compared to classical dynamic vegetation models. The authors just provide maps of height, NPP, AGB and AET that they compare visually with similar maps obtained from upscaled data or remote sensing to validate the model. This is not enough, especially that significant differences can be observed with the data. They should report statistical measures, such as correlation coefficients, RMSE and biases between modelled and observed maps.\n\nThe other point is that they are not validating model traits on present data (except height), which is the main outputs they use from the model in their study. They should compute at least at the level of a few sites (or regions) average traits values for LMA and deciduousness, from the genera or families present at these sites, using values from TRY, exactly as they do for the PETM.\n\nSome information should be added for several equations in the methodology section. For instance, equation (1) is not unit-independent. It is necessary to provide the units for SLA and all. Similarly, the units of the conductances should be provided in equation (9) (gm = 3.26 is a conductance and should have a unit). Also, in equations (18) to (22), the function N(alpha,0.05) should be clearly defined.\n\nIn equation (55), canopy conductance responds to (soil) water availability, but not to air relative humidity or saturation deficit, as in usual formulations of canopy conductance. This may need some comments.\n\nIt looks that the model neglects the diurnal cycle of photosynthesis and uses monthly average climatic data to calculate its rate. This is a strong approximation, that might be justified for computation efficiency. However, the impact on the absolute value of photosynthesis may be very important, because the process is not linear. Anyway, since relative differences are more important in this study, than absolute values, this hypothesis is probably fine, but it should be discussed.\n\nIn equation (57), Rn is daytime net radiation. It would be useful to provide some explanation on the way it is evaluated. It it just interpolated from the climate model outputs or does it react to the vegetation changes (e.g., through LAI impact on albedo) calculated by the vegetation model ?\n\nP 63, l 1282: km instead of km-1", "reviewer_3": "I thank the authors for this interesting and well-presented manuscript. The study is well designed, thorough and encompasses novel data compilation and modelling. The results are discussed coherently and extensively. Several key results stand out that are of high relevance to researchers in the fields of paleo and modern climate change and could have policy implications. To me, the key results are: a threshold warming beyond which terrestrial productivity is strongly reduced; vegetation changes during the PETM differed by latitude and productivity might have been stable at high latitudes; vegetation changes during the PETM would have contributed to the recorded d13C change across the event and could have controlled the strength and time scale of the silicate weathering feedback. These are all important findings that merit publication in a high profile journal. I suggest accepting it with minor revisions. It is a fascinating read and relevant to a wide audience.\n\nMy main comments are the following:\n\n1) I think this study contains multiple important results which could spelt out more explicitly and developed further:\nA - the threshold for exceeding dispersal strategy\nThe existence of a threshold warming for the capacity of vegetation to cope with warming is an important result in my view and of interest for the understanding of paleo events and the human-made climate change. How does this threshold depend on the speed of climate change? Rather than as an absolute temperature change, could this threshold be expressed as a warming rate?\nB - implications for the interpretation of the d13C anomalies during the PETM\nAs the authors mention multiple times, d13C changes in marine carbonates are a key feature of the geologic record of the PETM. It has even been used to estimate the magnitude of carbon fluxes during the event, usually without considering terrestrial d13C fluxes (e.g. Gutjahr et al. 2017, Komar and Zeebe, 2017, Haynes et al. 2020). The authors already allude to it but I think it would be interesting and could increase the relevance and reach of this paper by adding a brief discussion of the implications of the simulated terrestrial carbon stock changes in the present study would affect these estimates, and to put the simulated carbon fluxes in relation to those previously estimated.\n\n2) I think it would be a valuable addition to discuss how these results depend on the background state, e.g. how applicable they are to other climate events. The PETM starts from a much warmer background climate than today, with higher atmospheric CO2. How does this affect the potential of plants to adapt or migrate under warming compared to today? Does the CO2 fertilization effect tail off when atmospheric CO2 is above a certain threshold?\n\n3) Does the vegetation model consider disturbances like fire, storm or pest damage, which might have increased due to the PETM warming? If not, how might including these processes alter the results? Do they represent a best-case scenario?"}, "author_response": "We thank the reviewer for taking the time to review our manuscript, the positive feedback and the useful suggestions which have helped to improve the clarity of the manuscript and its conclusions.\n\nFirst, in a few places in the introduction (lines 46, 52, and 56) the text doesn’t clearly discriminate between biotic carbon stocks versus terrestrial carbon stocks specifically. Slight changes to the wording could clarify that existing d13C records do not uniquely implicate terrestrial organic carbon as a capacitor, as opposed to other organic carbon sources (marine organic carbon, hydrates, permafrost carbon, etc.). Changes in terrestrial organic carbon storage (vegetation and shallow soils, excluding permafrost) are one possibility, and this paper is evaluating the hypothesis that the terrestrial ecosystem response can account for the highlighted aspects of the global d13C records.\n\nWe thank the reviewer for ensuring clarity and we have revised the text to avoid implying that vegetation and soil carbon stocks are the only capacitor system of potential importance during the PETM. Particularly, we now mention other suggested capacitor systems (e.g., methane hydrates, permafrost) and clarified that we focus on the role of vegetation and soil carbon stocks (L43-50).\n\nNext, the Methods section answered most of the model-specific questions that I had while reading the text, but I think it would be helpful to include descriptions of the following in the main text (‘Trait based vegetation model’ section) in a third paragraph before Results:\n1) a brief description of how the vegetation model interacts with soils - e.g. the sentence in the Figure 1 caption that soil carbon is a function of net primary productivity and temperature-dependent heterotrophic respiration, 2) a brief description of how silicate weathering is calculated - e.g. move forward the description on Line 160 and clarification of whether soils and weathering are part of the vegetation model or computed offline, 3) a brief description of the experimental framework, e.g. how the time-dependent changes were constrained (50 kyr using pre-PETM constraints, 70 kyr of PETM constraints, then a linear CO2 decrease over 130 kyr with linear interpolation of climatic fields), and 4) a brief summary of the sensitivity tests - varying the adaptation and dispersal capacity and including CO2 fertilization and different degrees of warming. Specifically, how are different amounts of warming imposed? What is the amount of warming that is consistent with the CESM simulations (pre-PETM to PETM)? Is scaling applied to both spatial patterns in temperature change as well as the seasonal temperature range? These additional descriptions of both the model and the experimental framework would allow much more clarity in interpreting the subsequently presented results.\n\nWe have extended the methodological description in the main part of the manuscript as suggested by the reviewer (L103-145). We have clarified that soil carbon and silicate weathering are computed offline based on the results of the vegetation model and include details on their calculation in the main text. Further, we have added a third paragraph to the model description describing the climate boundary conditions used for the vegetation modelling, that the warming scenarios were derived by assuming different maximum CO2 concentrations during the PETM body, that the reference warming scenario from the original CESM simulations results in a 8°C land surface warming, and that monthly climate fields required by the vegetation model were derived by linear interpolation to the specified CO2 level from the available CESM simulations.\n\nA few more specific model questions (but that could go into the Methods):\nLine 99: are modern orbits assumed in calculating seasonality in both the pre-PETM and PETM simulations with CESM?\n\nThe CESM simulations used in this study were conducted considering a neutral orbital configuration with 0 eccentricity and 23.5° obliquity, which we have now clarified on L632-633.\n\nLine 523: what is the rationale behind selecting a 10 kyr model timestep?\n\nWe have clarified on L608-615 that the choice of the modelling time step in the vegetation model is arbitrary. Here, 10 kyr was chosen to have a reasonably abrupt onset of the PETM warming, while still being able to computationally cover the full 200 kyr length of the PETM and to conduct a large range of sensitivity tests regarding rate of evolution, dispersal, magnitude of warming and CO2 fertilization.\n\nLine 572: what is the temporal variation of climate input data? Monthly?\n\nWe have clarified in the main manuscript and the Methods section that the model uses monthly climate fields at each timestep (L135-136; L611-615)\n\nIn the discussion, it would be helpful to add some text to clarify the lack of feedback between the imposed climate boundary conditions and vegetation response and consider implications of that assumption. For instance, a modeled drop in biomass should generate not only a CO2 feedback via carbon release but also have a climate influence via changes in albedo, evapotranspiration, etc., and these cannot be represented with this experimental framework. Such feedback could potentially influence the timescale of the vegetation-mediated climate response.\n\nWe have clarified in the discussion (L389-399) that by imposing climatic boundary conditions, we do not consider such feedback with the climate system or the geological carbon cycle. We have further added to the discussion that the consideration of such interactions in a fully coupled manner will be important to understand how a potential productivity loss – as suggested by our results – may affect the magnitude and timescale of a PETM-like climate perturbation through positive feedback dynamics (e.g., release of biomass carbon).\n\nFurther, how have the results been shaped by assumptions in the modeling framework - specifically, the imposition of the PETM climate for 70 kyr? What if, instead, the experiment imposed a linear recovery in CO2 and temperature without a lag and more closely approximating the modeled response to a single pulse release of CO2 as has been simulated in models of the ocean-atmosphere (e.g. cGENIE, LOSCAR)? Would there still be a prolonged reduction in productivity, which in turn might explain why the actual PETM warmth persisted over a long duration? Given the number of sensitivity tests already included in the study - evaluating trait adaptation capacity, dispersal capacity, CO2 fertilization and degree of warming - I am not suggesting that it would necessarily add clarity to evaluate different assumptions about the climate boundary conditions, but it would be useful to discuss how the assumptions about these boundary conditions have impacted the results.\n\nWe thank the reviewer for highlighting this important discussion point, which we believe is key to communicate to the reader. Indeed, carbon cycle models like cGENIE and LOSCAR which focus strongly on ocean biogeochemistry generally model a fast onset of the climatic recovery after an initial pulse of carbon input. To obtain a prolonged warming as observed in proxy records (e.g., Fig 1A) they require additional carbon inputs during the PETM body (e.g., Penman et al. 2016, Zeebe et al. 2009). In these models, it is assumed that silicate weathering on land and subsequent carbonate sedimentation in the ocean effectively regulate climate as soon as temperatures rise. If we imposed a faster recovery of climate in our model, vegetation would recover faster as well, but there would still be a period of disrupted climate regulation for as long as the warming magnitude exceeds the adaptation capacity through dispersal and evolution. This temporary loss could then trigger positive feedback dynamics, including a loss of biospheric carbon storage and weathering, further prolonging the warming. We state in the discussion that the delay in weathering and organic carbon sequestration as a result of a vegetation disruption may offer an additional explanation for the delay in the climate regulation observed during the PETM (L365-385) and we now more clearly highlight (L389-399) that it may well be the response of vegetation that could determine the magnitude and duration of the PETM warming, as the reviewer suggests. To evaluate the effects of the vegetation disruption on the timescale of the climatic recovery through positive feedback dynamics, the proposed adaptation processes should be considered in a fully coupled manner with the climate system and the geologic carbon cycle. The main goal of the presented work and with the developed model was to evaluate the possibility that there was a widespread loss of vegetation functions during the PETM. Integrating these dynamics in a fully coupled Earth system and carbon cycle model over geologic timescales is an endeavour we hope to support with our findings, but that will require considerable further research and model development.\n\nFinally, in making comparisons between results for the PETM and the long-term future, do the authors have a sense whether the response is related to the relative or total amount of warming? Is 4 deg C a threshold because of an absolute temperature reached, or would 4 degrees of warming from a lower baseline have similar results, particularly because of the climatic niche constraint?\n\nWe have extended the discussion on the comparison between the PETM and anthropogenic climate change in line with the reviewer’s comments (L400-420). The vegetation response will likely depend on both the total and relative amount of warming. We now state that the absolute temperature limit before a loss of productivity occurs may lay higher for the present because we start in a colder background climate state. Consequently, the high latitude positive effects of warming, creating environments more optimal for photosynthesis (warmer, increased CO2) as we observed for the PETM, might extent into lower latitudes. However, the relative warming and the speed of warming are important, with current warming being considerably faster than during the PETM. This means less time for both dispersal and evolutionary adaptation. While both for the PETM and the present, slow evolution may induce an adaptive lag for a given relative climate shift, today, limited dispersal may have the potential to induce an additional adaptive lag because of the rapid speed of the climatic change.\n\nWe thank the reviewer for the positive feedback and for taking the time to go through the previous revision comments as well as for providing valuable feedback to ensure clarity of our methodological approach.\n\nFollowing the reviewer’s suggestion, we have split the previous Supplementary Figure 2 into four figures (Supplementary Figs. 2-5, referred to on L103), showing the visual comparison between modelled and observed values, a plot showing the spatial bias, and a scatterplot showing the modelled vs. observed values for height, NPP, AGB and AET, including the RMSE and correlation coefficient.\n\nWe have included an additional validation plot, Supplementary Figure 6 (referred to on L103), for a global LMA and deciduousness validation, where we compare the modelled latitudinal gradient in LMA for evergreen and deciduous plants compared to trait observations, as well as the current distribution of deciduous trees to the model predictions. Regarding LMA, we show that the model captures the main feature of an opposing LMA gradient for evergreens and deciduous trees with the cooling of climate towards higher latitudes as observed in field data. Regarding deciduousness, the model broadly covers the locations of Europe, Eastern North America and East Asia, where deciduousness is an abundant plant strategy (deciduous trees composing more than 50% of the vegetation composition).\n\nWe have revised Supplementary Table 1 and have included an additional Supplementary Table 2 to list all model parameters, constants and variables including their units (referred to on L451-452), and have clarified equations 18-22 (L572-575).\n\nWe thank the reviewer for pointing out this ambiguity in the methodological description. In the revised methodological description (L1195-1215) we have clarified that the canopy conductance does depend on both the local water availability as well as the atmospheric water demand. However, for the atmospheric water demand we follow the LPJ vegetation model (Schaphoff et al. 2018, Sitch et al. 2003), by applying Monteith’s empirical relation between evaporation efficiency and surface conductance and calculate the atmospheric demand using an equilibrium evapotranspiration rate as a function of temperature and radiation (Eq. 57). If the atmospheric demand exceeds the local water availability, canopy conductance and the canopy-atmosphere gas exchange parameter λ are recalculated in a bisection algorithm to find a conductance and photosynthetic rate that results in a closed water balance.\n\nWe have clarified that the vegetation model uses monthly average climate inputs in the main manuscript (L135-136). Further, we clarified in the Methods (L608-615) that the vegetation model is developed to obtain average vegetation dynamics across long, geologic timescales. To cover such timescales, the model considers multi-year monthly average climate inputs at every 10 kyr timestep for which it is run but does not cover any variability at a higher temporal resolution.\n\nIn the discussion, we have clarified that with the current model implementation there are no interactions between the vegetation model and the climate system (L391-395) and further clarified on L1209 that Rn is derived from the climate inputs.\n\nChanged accordingly, thank you (now L1338).\n\nWe thank the reviewer for taking the time to evaluate our manuscript, the positive feedback and the helpful comments to clearly communicate the implications of our findings to a broad audience.\n\nWe have extended the discussion on how the vegetation response depends on the absolute warming and how the difference in the rate of warming between the PETM and the present may affect the occurrence of an adaptive lag (L406-420). A direct comparison of warming rates between the PETM and current climate change is challenging, given the large uncertainty in the duration of the PETM onset, with estimates ranging from a few thousand years to 20 kyr (Turner, 2018). On the investigated warming timescales of 1-8°C over 10 kyr (now stated on L412-413), we show that dispersal was likely not a major limiting factor, and a species redistribution may have been able to compensate warming of up to 4°C land surface warming, after which climate adaptation through evolution became limiting. For the present, this absolute warming threshold may lay higher, as a colder background climate state likely results in the positive effects of warming and CO2 fertilization on productivity we observe for the PETM extending to lower latitudes. However, with carbon currently being released at a rate likely exceeding the PETM by around a magnitude, the speed of dispersal may become much more important at present for inducing an adaptive lag. In summary, we believe absolute warming caused the adaptive lag during the PETM, whereas at present, the rate of warming may be considerably more important, despite a potentially higher absolute warming threshold beyond which evolution becomes limiting.\n\nWe have extended and clarified the discussion on how our modelled terrestrial carbon stock changes compare to proxy-based carbon budgets and flux estimates for the PETM (L365-376). Recent studies suggest that approximately 5800 Pg of carbon were released during the PETM, that isotopically light carbon likely contributed substantially to the total carbon release (in addition to volcanic carbon and methane), and that a total organic carbon sequestration of 2000-2500 Pg of organic carbon would have been necessary to explain the δ13C records during the PETM recovery. Our simulations estimate a total carbon pool of 2000-2500 Pg of carbon in vegetation and soils and that 50% or more of this carbon could have been released during the PETM onset, followed by a lagged regrowth. Thus, according to our results, the terrestrial carbon dynamics may have contributed substantially to the global dynamics as derived from proxy records.\n\nWe have taken up the aspect of the background climate state in the comparison of the PETM and current climate change (L400-426). During the PETM, we show that increased temperatures and CO2 fertilization may have shifted conditions in the highest latitudes towards more optimal conditions for photosynthesis and productivity. In a generally colder climate and a lower initial atmospheric CO2 concentration, we expect these effects to extend further into lower latitudes. For a colder background climate state, we further expect that it requires more warming to have widespread regions were temperature and precipitation conditions are unsuited for primary productivity (e.g., temperatures exceeding critical thresholds for photosynthesis).\n\nWe thank the reviewer for highlighting this aspect. The current model does not consider such additional disturbances, which might have played a role during the PETM. We can expect that for a given rate of adaptation in the model, additional disturbances would lower the productivity and carbon sequestration and prolong the lag in effective climate regulation. We have added this aspect to the discussion (L339-340)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I reviewed a previous version of this manuscript. The authors have addressed all of my comments in both their response and revision. I think this is an excellent study and worthy of publication in its present form.", "reviewer_2": "Review of: “Loss of vegetation functions during the Paleocene-Eocene thermal maximum” by Rogger, J., Korasidis, V. A., Bowen, G. J., Shields, C. A., Gerya, T. V., and Pellissier, L., manuscript submitted to Nature Communications\nI read carefully the revised version of the manuscript. The authors addressed all the comments I made in my first review and modified the main manuscript and the supplementary text accordingly. From my point of view, the points raised by the other reviewers were also quite fully addressed. As a result, the paper can be published in its current form, after corrections of a few typos/grammatical errors I found while reading:\nPage 16, line 369: “a carbon cycle” instead of “a carbon carbon cycle”\nPage 17, line 411: “by around an order of magnitude” instead of “by around a magnitude”\nPage 20, line 462: “are used to describe” instead of “are used describe”\nPage 48, Legend of Supplementary Figure 6, A & B: the figures show LMA versus latitude, while the legend refers to SLA. All occurrences of SLA in the legend should be replaced by LMA. Also, legend for figure D is missing.", "reviewer_3": "I thank the author team for their efforts in revising the manuscript.\nI re-read the revised manuscript a few times now and think that, while future modelling studies should take more of the complexity of the coupled atmosphere-land-ocean system into account, this current manuscript presents a well rounded study and will surely spark important debates. I look forward to seeing this manuscript published."}, "author_response": "We thank the reviewer for the helpful comments during the review process.\n\nChanged accordingly.\n\nChanged accordingly.\n\nChanged accordingly.\n\nThe legend of Supplementary Figure 6 has been revised accordingly, and a description of panel D has been added.\n\nWe thank the reviewer for the helpful comments during the review process.\n\nWe thank the reviewer for the helpful comments during the review process.\n\nWe have checked the link, and all code files and data are accessible under:\nhttps://zenodo.org/records/17252973"}]

Decision letter