Low atmospheric CO2 levels before the rise of forested ecosystems – Watts Up With That?

Nature Communications volume 13, Article number: 7616 (2022) Cite this article

Abstract

The emergence of forests on Earth (~385 million years ago, Ma)¹ has been linked to an order-of-magnitude decline in atmospheric CO₂ levels and global climatic cooling by altering continental weathering processes, but observational constraints on atmospheric CO₂ before the rise of forests carry large, often unbound, uncertainties. Here, we calibrate a mechanistic model for gas exchange in modern lycophytes and constrain atmospheric CO₂ levels 410–380 Ma from related fossilized plants with bound uncertainties of approximately ±100 ppm (1 sd). We find that the atmosphere contained ~525–715 ppm CO₂ before continents were afforested, and that Earth was partially glaciated according to a palaeoclimate model. A process-driven biogeochemical model (COPSE) shows the appearance of trees with deep roots did not dramatically enhance atmospheric CO₂ removal. Rather, shallow-rooted vascular ecosystems could have simultaneously caused abrupt atmospheric oxygenation and climatic cooling long before the rise of forests, although earlier CO₂ levels are still unknown.

Introduction

Atmospheric CO₂ is a greenhouse gas that has affected Earth’s climate throughout geological history^2,3. Its variation in the past informs us about the natural long-term sources and sinks. In the absence of anthropogenic fossil fuel combustion, the dominant atmospheric CO₂ source is volcanic outgassing, and this source is balanced mainly by the removal that occurs when CO₂-bearing fluids chemically react and weather silicate rocks followed by deposition of carbonate in the oceans⁴. The dissolution of silicate minerals in the weathering zone occurs via interactions between the terrestrial ecosystem and geological processes that make fresh rock available at the surface for reaction. Yet, the role of biology and the CO₂-sensitivity of the feedbacks governing global CO₂ removal is debated^5,6,7,8,9. Enhanced continental weathering is suggested to have caused a decline in atmospheric CO₂ pressure (pCO₂) from a level ~10 times higher than today’s concentration^3,10,11 to near modern levels linked to the Devonian-Carboniferous transition from greenhouse to icehouse conditions in response to the afforestation of the continents. This process is accompanied by burial and preservation of organic matter that also influence atmospheric CO₂ levels and, in turn, acts as the main long-term source of atmospheric O₂. However, recent geochemical evidence and Earth system models¹² suggests atmospheric oxygenation occurred well before trees evolved on the continents ~393–383 Ma^1,3. Further, the temporal correlation between plant colonization and the Permo-Carboniferous glaciation has been disputed⁸. There is compelling evidence that Earth also transitioned into a glaciated state in the Ordovician-Silurian¹³. Yet, the link between glaciation and atmospheric CO₂ is complicated, and palaeoclimate models shows that glaciations could persist even at 12–14 times pre-industrial atmospheric levels (PIAL, 280 ppmv)¹⁴. Therefore, a precise reconstruction of atmospheric pCO₂ in relation to plant evolution is key to assess the impact of the terrestrial biota on Earth’s climate. Here, we show that atmospheric pCO₂ was markedly lower than previously thought when trees and forests appeared on our planet.

In the canonical view, atmospheric CO₂ concentrations were one order of magnitude above pre-industrial levels in the early Palaeozoic^{3,10,15,16,17}, although more recent studies suggest levels much closer to today¹⁸. Traditionally, palaeo-CO₂ estimates from proxy data come with large and sometimes unbounded uncertainty^19,20. For example, in the Late Ordovician (~445 Ma), there is evidence for high pCO₂ levels of 17 ± 4 PIAL (1 sd, standard deviation) from CO₂ hosted in pedogenic goethite from the Neda Formation in Wisconsin, USA¹⁹. This paleosol shows a coupling between CO₂ content in goethite and its isotope composition interpreted to reflect variable mixing of atmospheric CO₂ and soil respired carbon substituted into the goethite mineral lattice. The range above reflects all analytical errors propagated through the calculation (see details in the supplementary information, SI), and we note that the error could still be larger if the CO₂ surface adsorption properties²¹ on natural goethite deviate from that of phosphated goethite grown in the laboratory²². This has never been verified in modern soils, and the proxy has also never been applied at any other time in Earth history.

Further, a systematic decline of atmospheric CO₂ levels through the Devonian from 5 to 0.7 PIAL has been inferred from the carbon isotope compositions of pedogenic carbonate^2,16,17. Although the declining trend may be real, the absolute atmospheric CO₂ level reported from pedogenic carbonates from this time interval have been adjusted further down²³ and are likely systematically overestimated because of lower productivity in early Palaeozoic soils relative to modern soils¹⁷. Propagating the uncertainty associated with assumed model parameters that cannot be independently constrained from the rock record, such as the proportion of soil-respired CO₂ in the soil and its isotope signature, shows that the absolute palaeo-pCO₂ estimate obtained this way carries one order of magnitude uncertainty (see supplementary Fig. 15 and the supplementary information for details).

In addition, the first reported evidence for high Palaeozoic CO₂ levels (~16 PIAL)²⁰ comes from low stomatal density (mm⁻²) in some fossil plants (i.e. Aglaophyton, Sawdonia). These anatomical features were interpreted as evidence that these plants had adapted to minimize water loss in a high CO₂ atmosphere, but it is essential to compare to plants with similar gas exchange anatomy and behavioural control over water loss rate and CO₂ uptake. This becomes very problematic with this group of enigmatic early vascular plants that lack living descendants. In comparison fossils of lycophytes that co-occur with these extinct plant groups that do have modern relatives with similar physiology display similar stomatal density as their modern descendants (Supplementary table 5) suggesting that the high CO₂ predictions inferred from Aglaophyton and Sawdonia specimens could be erroneous and that early Palaeozoic atmospheric CO₂ was much closer to the modern level²⁴.

Recently, atmospheric CO₂ levels were found to be only modestly elevated in the late Ordovician (~400–700 ppm) and Mid-Devonian (~700–1400 ppm) based on carbon isotope data of marine phytoplankton recorded via a diagenetic product of chlorophyll, phytane¹⁸. No phytane-CO₂ data is reported between ~432 Ma and 390 Ma (supplementary Figs. 19–21). The phytane proxy shows the expected response at elevated CO₂ concentrations and predicts declining CO₂ levels as a function of distance to a modern CO₂ seep²⁵ (albeit with larger errors than reported from Palaeozoic phytane). When applied to modern phytoplankton, the phytane proxy predicts a wide range of atmospheric CO₂ levels even today (~300–1200 ppm)²⁶. The accuracy of this proxy depends on variables that are not easily detectable from the geological record, including the ratio of atmospheric CO₂ to dissolved CO₂ of the seawater in which the phytoplankton grew, the nature of the phytoplankton species, other sources of phytane, the growth rate of phytoplankton and, thus, the isotopic response to ambient dissolved CO₂ concentration of the phytoplankton from which phytane was derived²⁶. Given the many uncertainties, the errors associated with the phytane CO₂-proxy are likely large.

By the same fundamental principle plus an opportunity to constrain all involved parameters, carbon isotope fractionation in plant tissue from terrestrial plants (Δ_leaf) is sensitive to the internal CO₂ concentration in the substomatal cavities that, in turn, depends on water availability and ambient pCO₂ in the environment²⁷. Recently, the combination of isotope data and stomatal parameters has led to the development of a mechanistic proxy for prediction of palaeo-pCO₂ based around leaf-gas exchange²⁷. Applying this approach to the fossil record shows that post-Devonian atmospheric CO₂ levels were <1000 ppm most of the time, but CO₂ estimates from the Lower and Middle Devonian using this model show considerable variation with estimates ranging from 530–2853 ppm with one outlier reaching 7320 ppm²⁷.

Importantly, this mechanistic proxy has also allowed for both the accuracy and the precision of the prediction to be tested allowing for the development of realistic constraints around the predicted palaeo-pCO₂ to be assessed. Thus, the framework that underpins the proxy allows for the re-evaluation of CO₂ predictions from both long-term carbon cycle models which have large uncertainty bracketing their predictions^11,12,15 and predictions from first generation palaeo-pCO₂ proxies that generated estimates in pCO₂, which were largely unconstrained.

The mechanistic model by Franks et al.²⁷ uses both carbon isotope data and stomatal features (density and size) obtained from the fossils to predict palaeo-pCO₂. Both characteristics are an expression of how the plant, when alive, was adapted to its local environment. Consequently, there is a strong relationship between these different “recorders” of plant climate interactions. For example, the carbon isotope signature of leaf tissue has been demonstrated to be strongly influenced by water availability²⁸ and the utility of mechanistic proxies underpinned solely by isotope data^29,30,31 to be questioned^{32,33,34,35,36}. Due to the co-dependency of isotope and stomatal parameters it is recommended that these data are ideally extracted from the same sample whenever possible or from contemporaneous sedimentary deposits. In taking this approach, the accuracy of the prediction should be improved because it simultaneously accounts for variable CO₂ assimilation rate and water conductance in response to changing water availability²⁷. Here, we calibrate a leaf gas-exchange model for lycophytes and estimate atmospheric CO₂ levels with uncertainties by applying it to some of the oldest representatives found in the fossil record.

Results and discussion

Calibration of a robust palaeo-CO₂ barometer for lycophytes

To anchor the early atmospheric CO₂ records, we calibrated a mechanistic leaf-gas exchange model²⁷ in living representatives of the most ancient vascular plant taxa (Lycopodiaceae) constrained by their stomatal density, stomata pore size, and carbon isotope compositions, and applied it to fossil remains of Lower and Middle Devonian lycophytes with the same stomatal anatomy³⁷. When trying to use the plant fossil record to estimate lower Palaeozoic CO₂ levels at the time when vascular plants underwent their initial radiation, it becomes challenging to use paired sets of stomatal and isotope data. This situation arises due to the mode of fossil preservation, the need for destructive analysis for carbon isotopes, and the intrinsic scientific value of curated specimens. Therefore, we have further developed the mechanistic model²⁷ to incorporate a solver routine that allows for the propagation of a wide range for each of the observable parameters required to deliver palaeo-pCO₂ estimates. This approach has been developed as it allows for all variation between local environments and/or a temporal mismatch between data sources to be folded into the model predictions. Again, to ground truth that our solver approach delivers reasonable estimates of palaeo-pCO₂, we have benchmarked the approach by applying the technique to our extant lycophyte dataset.

Initially, we calibrated the gas exchange model²⁷ using two lycophyte species (Huperzia phlegmaria and H. squarrosa) presumed physiologically similar to the Devonian lycophytes and grown at known ambient CO₂ levels for ~8 years in natural light within a soil substrate and under optimal Relative Humidity (RH) of ~80% in a greenhouse of the Botanical Garden in Copenhagen. Measurements of the leaf carbon isotope fractionation (Δ_leaf), stomatal density (SD), and stomatal pore length (p) (Supplementary table 2; Supplementary Data 1) suggest that their CO₂ assimilation rate at reference CO₂ (A₀) was indistinguishable from that of other modern lycophytes (A₀ = 3.7 ± 1.6 μmol m⁻² s⁻¹, 6 species). Using this average A₀ value as representative for lycophytes, and other parameters fixed, such as the operational stomatal conductance efficiency (ζ) with small effect on pCO₂ estimates (see supplementary text for details), we were able to reproduce measured CO₂ levels of 448 ± 51 ppm (daily variation, 1 sd) in the greenhouse with a CO₂ prediction of 407+24−22407−22+24 ppm and 545+69−56545−56+69 ppm (median ± quartiles) for the two distinct species of lycophytes, respectively.

Following this initial screening we further explored the mechanistic model to derive estimates of CO₂ by developing a dataset based on non-paired data; i.e., where stomata and isotope data comes from H. phlegmaria and H. squarrosa from different growth conditions to mimic fossil samples that come from geographically/environmentally distinct sedimentary deposits of similar age. Based on cultivated lycophytes grown under dry and humid conditions, we find that unpaired data could carry up to ~100 ppm additional uncertainty, if stomatal and isotopic data were derived from plants that lived under vastly different conditions; i.e., stomatal data from slow-growing or even semi-epiphytic plants under extremely dry conditions with isotope data from plants grown under optimal, humid conditions or vice versa (Supplementary Data 1). When explored using this data framework, the revised predictions of pCO₂ shows good agreement with both the actual growth CO₂ and values of CO₂ predicted using the standard model parameterization (Figs. 1–2; Supplementary Data 1). Thus, the simultaneous analyses of fossil lycophyte Δ_leaf, SD, and p should deliver both accurate and precise predictions of palaeo-pCO₂ irrespective of whether the data are generated from paired samples (stomatal and isotope data from the same locality) or from different localities. Consequently, by applying the mechanistic model²⁷ to lycophyte data, we are now able for the first time to resolve atmospheric pCO₂ accurately and precisely prior to the emergence of forests.

Fig. 1: Atmospheric pCO₂ levels are derived from leaf carbon isotope fractionation (Δ_leaf), Stomata Density (SD), and stomata pore length (p).

a Δ_leaf is determined from the carbon isotope composition (δ¹³C) of the fossil flora (green), marine carbonates (blue) and calculated atmospheric CO₂ (light blue) (Supplementary Table 4). b Plants grown under higher ambient CO₂ levels yields a higher Δ_leaf and/or lower stomata density. The curves represent gas-exchange model calibrated for modern lycophytes with a reference CO₂ assimilation rate at modern atmospheric levels (A₀) of 3.5 µmol m⁻² s⁻¹, an operational stomatal conductance efficiency (ratio of operational to maximal stomatal conductance) ζ = 0.2, stomata pore shape β = 0.6 and p = 18 ± 2 µm comparable to ancient and modern lycophytes (see supplement for details). Here, a fossil lycophyte from the Devonian is used as an example. Devonian lycophytes typically display Δ_leaf values close to 20‰ and stomata densities between 15 and 25 mm⁻¹ (Supplementary Data 1).

Fig. 2: Error propagation of the mechanistic pCO₂ proxy applied to Givetian Drepanophycus sp. from the Hamilton Group, NY, USA.

Atmospheric CO₂ levels at 410.8–382.7 Ma (Pragian through Givetian) were reconstructed from 66 fossil lycophytes representing three distinct genera (Asteroxylon, Baragwanathia and Drepanophycus) from 13 geological deposits at nine geographically distinct localities. From this, we derive four pCO₂ estimates from paired data and six predictions from unpaired data (Supplementary Data 1 for location and data details). However, we note that the stomata densities from all localities were remarkably similar (±1.5 mm⁻², 1 sd) and there is a degree of similarity in Δ_leaf (18.3 ± 1.1‰, 1 sd). Using paired data pCO₂ predictions are constrained to between 525+139−101525−101+139 ppm and 695+99−73695−73+99 ppm (median ± quartiles of the probability distribution). Unpaired predictions give a range of pCO₂ estimates constrained to between 532+77−76532−76+77 ppm and 715+140−102715−102+140 ppm. Our data yields consistent estimates of CO₂ regardless of the method used to compile the raw data supporting the predictions (Fig. 3; Supplementary Data 1). Thus, we conclude that these data substantiate that atmospheric CO₂ levels were only 1.9–2.6 times above pre-industrial levels during the ~30-million-year time interval when plants evolved tree stature and forests appeared on Earth.

Fig. 3: Summary of plant evolution and evolving atmospheric composition versus age.

Atmospheric pCO₂ constraints from 410–380 Ma lycophytes (Asteroxylon, Baragwanathia, and Drepanophycus) with errorbars smaller than the size of the circles. Paired stomata and isotope data from the same locality are highlighted with black circles. Atmospheric pO₂ is constrained by charcoal evidence for wildfire since ~423 Ma54, 55 (red horizontal line) and fossil roots (red circles)⁵⁶ that sets minimum pO₂ levels according to controlled calibrations in the laboratory^56, 84. Modelled evolution of atmospheric pCO₂ (blue curve) and atmospheric pO₂ (red curve) derived by adjusting the COPSE Reloaded model for the coupled biogeochemical cycles such that continental weathering efficiency of early vegetated ecosystems scales with the physical weathering record of fines in terrestrial deposits⁴⁶ and outgassing from the Earth scales with continental arc volcanism. The ‘COPSE reloaded’ model are shown for comparison (thin, dashed curves)¹². For further details, see supplementary information. The emergence of land plants, vascular plants, arborescent plants with deep root systems and seed plants are shown with thin lines representing their origin by molecular clock estimates and thicker arrows representing fossil occurrences⁴⁵. Llan Llandovery, W Wenlock, L Ludlow, P Pridoli. The uncertainty of absolute age assignments are defined biostratigraphically, approximately ±1.5 Ma.

Two potentially confounding factors on the magnitude of leaf carbon isotope fractionation are humidity and atmospheric O₂ levels. First, plants down-regulate stomatal conductance and increase their water utilization efficiency in drier habitats38. Lycophytes have passive stomata control and can adapt in drier habitats by minimizing stomata density and size and/or survive at lower CO₂ assimilation rate, which would then be expressed in lower Δ_leaf²⁷. We verified this adaptation experimentally using H. squarrosa that grew semi-epiphytically (Supplementary Fig. 6) under drier conditions (RH~60%) resulting in significantly lower Δ_leaf (13.3 ± 0.2‰) than the Δ_leaf (20.1 ± 0.9‰) of plants grown under optimal humidity (RH~80%; Supplementary table 2). This effect is also observed on the isotope composition of natural populations of C3 plants, but is significant only in areas where the mean annual rainfall is below ~1000 mm/yr (Supplementary Fig. 10)²⁸. In the drier greenhouse, H. squarrosa plants had a lower, albeit more variable, stomata density (16.1 ± 4.7 vs. 20.1 ± 2.8 mm⁻²) and a similar pore length (22.8 ± 1.0 µm vs. 21.6 ± 2.7 µm). Although, we do not know that these lycophytes had fully adapted to the drier greenhouse conditions, the pCO₂ estimate derived from such plant material yields an under prediction of ambient glasshouse pCO₂ level by -128+79−67128−67+79 ppm, which could be explained by a slightly lower operational to maximal stomatal conductance ratio (ζ) than plants grown under ideal (natural habitat) conditions (e.g. ζ of 0.14 vs. 0.20; see supplementary information for details).

To assess this humidity effect further, we used an atmosphere-ocean palaeoclimate model of intermediate complexity (CLIMBER-3α) to evaluate the Devonian climate with 500 ppm of atmospheric CO₂, at the lower end of the range of our new CO₂ estimate. The model was run to better constrain conditions where the early lycophytes lived (see supplementary information for details) and shows a temperate planet with mean tropical surface air temperatures of 24.1–24.6 °C. The lycophyte floras in Australia, Germany, and China grew in the monsoonal belt with high mean annual rainfall (MAR) exceeding ~1000 mm/yr. But, Asteroxylon (Rhynie Chert, Scotland) and Drepanophycus floras from Maine, USA; New Brunswick and Gaspé, Canada are found in the slightly drier subtropical zones. Due to its coarse spatial resolution, our model cannot accurately capture local variations in rainfall and humidity at the fossil sites, but the model outputs do suggest growth in either humid or relatively humid environments. Importantly, we see no coupling between Δ_leaf or predicted pCO₂ from each fossil localities and the modelled MAR, RH%, or palaeolatitude. Also, the Δ_leaf data is not suggestive of arid growth conditions. Given these data, there is no indication that the Devonian lycophytes lived under semi-arid/ arid conditions or the palaeo-pCO₂ proxy reported here should deliver a biased prediction as a function of a down regulation in stomatal conductance due to growth in an arid environment.

Secondly, atmospheric O₂ can also affect the leaf carbon isotope fractionation because photosynthetic CO₂ fixation competes with photosynthetic O₂ fixation on the Rubisco enzyme³⁹. This effect is more pronounced in modern plants grown at sub-ambient CO₂ and super-ambient O₂ levels⁴⁰. The Devonian atmospheric pO₂ levels were likely ~15 to 20 atm%¹² and not much lower than today, so any O₂ effect on our palaeo-pCO₂ estimates will likely be small. Controlled growth experiments with vascular plants, including one species of lycophytes (Selaginella kraussiana), show a small positive shift in Δ_leaf (0.5 ± 0.4‰) when grown at sub-ambient O₂ levels (16 atm%)⁴⁰. This effect has been ascribed to O₂-sensitivity of the photorespiratory compensation point⁴¹ and has a negligible effect on past atmospheric CO₂ estimates compared to the O₂-insensitive parameterisation of the mechanistic proxy (Supplementary table 4).

Importantly, our new and more precise CO₂ predictions are radically lower than previous reported mean values, but they are within uncertainty envelope of other pCO₂ proxies, including the pedogenic carbonate record when considering that Devonian soils were arguably less productive with lower soil respiration rates¹⁷ and lower CO₂ concentration than modern soils (supplementary Fig. 14).

Reconfiguring the leaf-gas exchange model to solve for CO₂ assimilation rate (A₀) allows us to explore ecophysiological performance of early land plants that have no living relatives. This approach enables us to contextualize the very low stomatal density of these plants. Model results suggest that these plants had substantially lower CO₂ assimilation rates than contemporaneous lycophytes (A₀ < 3.7 μmol m⁻² s⁻¹); Aglaophyton (0.40 ± 0.10 μmol m⁻² s⁻¹), Rhynia (0.62 ± 0.13 μmol m⁻² s⁻¹), Horneophyton (0.84 ± 0.21 μmol m⁻² s⁻¹), Sawdonia (0.89 ± 0.04 μmol m⁻² s⁻¹) and Nothia (1.27 ± 0.29 μmol m⁻² s⁻¹) (Supplementary Data 1). Thus, it is plausible to suggest that the low rates of assimilation in this grouping of plants could be a factor in their eventual displacement. The calculation of low assimilation rates suggest that the sporophyte could be physiologically tied to the gametophyte, suggestive of a degree of matrotrophy, as suggested for some species of Cooksonia⁴².

Atmospheric CO₂ and early afforestation

The atmospheric pCO₂ curve (Fig. 3) illuminates how the early history of land plants and their colonization of the Earth’s land surface affected the Earth’s climate system. The origin of trees in the Mid-Devonian (393-383 Ma)¹ represents a major evolutionary shift in which arborescence evolved independently in three distinct plant lineages (lycopsids, sphenopsids, ferns). Progymnosperm trees towering up to ~30 m diversified and became widespread from palaeo-equatorial to palaeo-boreal latitudes and from seasonally dry to wet habitats^1,43. In contrast to the earlier lycophyte-dominated flora of typically ~0.1 m in height and with rhizoids penetrating only few cm into the ground, the evolving trees had both deeper roots for anchoring and a highly advanced root system akin to modern seed plants¹. Still, atmospheric CO₂ levels remained rather low and may have declined¹⁷ by, at most, a few hundred ppm during the Devonian emergence of forests.

The absolute magnitude of the early Palaeozoic pCO₂ decline is still poorly constrained, and there is a conflict between some low pCO₂ estimates obtained from Mid-Ordovician marine phytane records and the high pCO₂ estimates from pedogenic goethite¹⁸. If we trust that the Late Ordovician atmosphere had a higher pCO₂ than today, at least outside glacial maxima (supplementary Fig. 12), our new results points to a dramatic and abrupt decline from ~9 (unbound) PIAL to ~1.9 ± 0.3 PIAL that took place within a relatively short time interval ~445–410 Ma before forests appeared. Coincident with this, shallow shrub-like vascular plants spread on the continents (e.g. the Eophytic flora⁴⁴) and a dramatic shift in the physical weathering regime driven by the evolving terrestrial ecosystems is recorded by the retention of fine grained sediment in continental deposits^45,46,47. We propose that the earliest vascular vegetation promoted the exposure of more mineral surface area to weathering fluids and amplified global silicate weathering on the continents (far more than subsequent deep-rooted ecosystems could do) owing to a greater nutrient loss from less developed soils, and therefore, forcing a higher weathering demand.

To simulate the effect of the establishment of early vascular plant ecosystems, we used a dynamic model (COPSE Reloaded) for the coupled biogeochemical cycles to predict atmospheric pCO₂ and pO₂ trajectories (see Methods). Our model predicts a massive atmospheric CO₂ decline from ~2500 ppm to ~500 ppm in only ~30 Myr in response to enhanced silicate weathering by early vascular plants (Fig. 3). To exemplify such a scenario, we updated several forcing functions in the latest COPSE model^12,48,49 (section S5; Supplementary Figs. 19–21). The weathering forcing (W) was adjusted to scale in proportion to the plant-induced effect on mudrock retention in continental deposits normalized to the Carboniferous average⁴⁶. This is justified because mineral surface area is a key factor facilitating mineral dissolution during chemical weathering. Also, we scaled up plant evolution in concert with the radiation of vascular plants rather than non-vascular plants in previous COPSE models in an attempt to both capture a more extensive plant coverage^50,51 and also accounting for selective P weathering to mimic a greater weathering demand of early vascular plants with primitive root-like systems^49,50,51. Lastly, volcanic outgassing rates in the Lower and Middle Devonian were adjusted so that the CO₂ flux emanating from Earth’s interior scales with subducted carbonate platforms rather than to global seafloor spreading rates as was assumed in previous models^11,52. This was done by assuming a linear relationship between outgassing rate, continental arc volcanism and the relative abundance of young to older arc-derived grains in sedimentary deposits⁵³. The Late Silurian to Middle Devonian (~430–370 Ma) outgassing rates is then ~1.0–1.2 times modern levels (in contrast to ~1.5 in previous models) characteristic of ages when supercontinents assembly.

The revised model simultaneously predicts that the evolving vascular ecosystems also caused a rise in atmospheric pO₂ consistent with multiple lines of geochemical evidence for Earth’s oxygenation¹² including wild fire evidence supporting atmospheric O₂ levels above 15 atm% since ~420 Ma^54,55 and fossil roots⁵⁶ (Fig. 3). Previous models have suggested a two-stage transition with relatively high pCO₂ and high pO₂ in the Silurian¹², but our new data and model offers a simpler solution where a synchronous atmospheric O₂ rise and CO₂ decline to near-modern conditions happened in the same event.

Ultimately, the composition of Earth’s atmosphere is governed by an interplay between biological and geological processes and how land plants and their root symbionts affect the physical and chemical weathering processes on land^45,57. The difference in atmospheric CO₂ by forestation post-vascular colonization was at most a few hundred ppm⁵⁸. Thus, global climatic cooling is not strongly linked to afforestation⁵⁹, but rather to how the terrestrial biosphere extracts and maintains nutrients from its planetary substrate.

The Mid-Palaeozoic climate revisited

The Devonian climate has previously been described mostly as a warm greenhouse that transited into a colder state with polar glaciations in the Late Devonian-Carboniferous^10,15. At only 500 ppm CO₂, however, our Early Devonian palaeoclimate model predicts a clement climate with global mean annual temperatures of 14.8–15.7 °C for different orbital configurations and a strong latitudinal temperature gradient not too dissimilar from today. Mean tropical surface air temperatures are 24.1–24.6 °C. The predicted sea surface temperatures in tropical Gondwana, South China, Western Laurentia, and Europe are consistent with Early and Middle Devonian palaeotemperature proxy records based on oxygen isotope compositions of phosphatic conodonts and well-preserved calcitic brachiopods (Supplementary Data 2)^60,61. Our paleoclimate model with 500 ppm atmospheric CO₂ predicts significantly lower temperatures in the subtropical and temperature zone where there are currently no precise palaeotemperature estimates (see supplementary information section S4.2 for further discussion)⁶². Polar sea ice and snowfall on Gondwana are predicted during winters (Fig. 4, Supplementary Figs. 14–16). CLIMBER-3α is not coupled to a dynamic ice sheet model, so the extent and persistence of ice sheets is indirectly constrained by the extent of snow and sea ice. We find that the Earth system was climatically stable and that a cascade into a snowball-style glaciation due to the ice-albedo feedback is highly unlikely at this low CO₂ level even under the weaker solar insolation, a finding that is also supported by model simulations bracketing the Early Devonian^63,64,65. Therefore, our results suggest that Earth’s climate was in an icehouse state with partial glaciation on the south polar continent Gondwana above 60-80°S palaeolatitude (Fig. 4). We suggest from modern palaeoclimate models^14,66,67 that the original conjecture¹⁰ that several thousands of ppm CO₂ in the atmosphere were necessary to compensate for the ~3% weaker Palaeozoic Sun should be abandoned, in favour of much higher climate sensitivity of CO₂ as supported by the new palaeo-CO₂ proxy data³.

Fig. 4: Early Devonian (415 Ma) annual mean surface-air temperature for an atmospheric pCO₂ of 500 ppmv, modelled with the coupled climate model CLIMBER-3α.

The palaeogeographical configuration of Scotese at 415 Ma83 is used, pre-industrial levels of other greenhouse gasses and idealized orbital parameters are assumed (circular orbit, obliquity of 23.5°; see supplementary Fig. 17 for two different orbital states). See supplement for details about modelled precipitation and humidity. Nine localities are marked where the studied fossil lycophyte flora have been found in 13 distinct outcrops: 1 – Green County, NY, USA and Schoharia Co, NY, USA; 2 – Traveller Mountain, ME, USA; 3 – North Shore and Seal Rock, Gaspe Qbc., Canada; 4 – Abitibi River, Ontario, Canada; 5 – Dalhousie, New Brunswick, Canada (close to Maple Green); 6 – Munchshecke, Siegburg and Mosel Valley, Germany; 7 – Rhynie Chert, Aberdeenshire, Scotland; 8– Yea district, Victoria, Australia; 9 – Xinjan, China. Six localities with palaeo-SST data from calcitic brachiopods are marked with white stars: a – Section Madène el Mrakib, Anti-Atlas Mountains, Morocco; b – Colle area, Cantabrian Mountains, Spain; c – Eifel Mountains, Germany; d – Buffalo Quarry, Iowa, USA; e– Mafeking Quarry, Manitoba, Canada; f– Gur’evsk town, Altai Mts./Salair area, Siberia. Seven localities with palaeo-SST data from phosphatic conodonts are marked with purple stars: g – Changputang, SE Yunnan, China and Nayi, Guangxi, China; h – Broken River, Queensland, Australia; i – Buchan, Boola Quarry, Tyers Quarry, Wellington, Victoria, Australia; j – Mungallala and Windellama, Wellington New South Wales, Australia; k – Barrandov and Na Skrabku, Prague basin, Czech Republic; m – Nevada, USA; n – CaiZiyan, Guangxi, China. Phosphatic SST constraints from d –Glory Quarry, Iowa, USA and b – Puech de la Suque, France are plotted along with the calcitic SST data from nearby outcrops.

Direct evidence for glaciation occurs in Late Devonian – Carboniferous deposits68, but is rare in Lower- and Middle Devonian strata. Indirect evidence of eustatic sea level change occurs in the Lower Devonian sedimentary successions from North Africa that could potentially have a glaciogenic origin⁶⁹. Further, oxygen isotopic palaeotemperature records⁶⁰ from the shallow tropical and subtropical oceans are compatible with our predicted palaeoclimate data, but also display considerable local temperature variations⁶¹. Also, Devonian plants are mostly found at lower-mid palaeolatitudes (~45°S), where the continent is predicted permanently ice-free. One exception might be the early Lockhovian Cooksonia-rich flora from the Paraná basin in Brazil positioned at high palaeolatitudes (~70°S)⁷⁰. However, even at 500 ppm CO₂, we find that the snow cover on Gondwana was not always perennial (supplementary Fig. 15). Although, the Paraná flora notably predates our record (~419 Ma) when atmospheric CO₂ levels were presumably higher, plants could also have grown at high latitudes either during warmer Southern summers or because climate was warmer for other reasons than high CO₂ (e.g. higher CH₄, N₂O).

In summary, the revised atmospheric CO₂ record is compatible with all palaeoclimatic constraints from the geological record, with global palaeoclimate models, and with global models for the coupled biogeochemical cycles that include a larger weathering demand of the earliest terrestrial flora and where early afforestation played only a minor role on global CO₂ sequestration.

Methods

Samples

Carbon isotopes, stomata density and stomata pore size from fossil leaves was compiled from the literature (see supplementary information) along with new data collected from the Baragwanathia flora in Victoria, Australia. Plant macrofossils from the fossil collection of Museums Victoria were originally collected by Isabel Cookson from Mt Pleasant Rd, Victoria, Australia⁷¹. The fossils are preserved as incrustations in fine-grained sandstone sandwiched in a 130 m thick stratigraphic section of mostly siltstone and shale (supplementary Fig. 1), interpreted as occasional bursts of high-energy turbidites carrying allochthonous fossils from shallower waters into very-low-energy marine depositional environment⁷². Index fossils (Uncinatograptus sp. cf. U. thomasi and Nowakia sp. ex gr. N. acuaria) confines the flora to the Pragian or earliest Emsian, corresponding to ~409.1 ± 1.5 Ma according to GTS2020⁷³.

Small fossil fragments found in four specimens (#15154, #15173, #15174, #15183) which also contain larger fragments of Baragwanathia longifolia and Zosterophyllum australianum, were selected for analyses based on the preservation of black organic matter contained within the brown mineralized fossils (supplementary Figs. 2–4). From each fragment (37 in total), 0.05–7.34 mg of material was extracted using a scalpel or a 0.8 mm Dremel drill. Only material visibly containing black organic matter was extracted.

Carbon isotope analyses

Carbon isotopic analyses for were performed in the Syracuse University GAPP Lab using an automated ‘nano-EA’ system adapted from that described in Polissar et al.⁷⁴. The Syracuse University nano-EA comprises an Elementar Isotope Cube elemental analyzer coupled to an Isoprime 100 continuous-flow stable isotope mass spectrometer via an Isoprime Trace Gas analyzer. Though the presence of carbonate was not suspected, nor observed via testing of sample powders with 6 N HCl under a binocular scope, we decided to decarbonate the sample materials prior to analysis to ensure that sample materials were carbonate free. Sample powders were fumigated in ashed glass vials in the presence of neat hydrochloric acid within an evacuated glass bell jar for 24 hours and then were dried in an oven at 40 °C. For isotopic analysis, sample materials were transferred to silver or quartz cups (6 ×6 mm; EA Consumables) and nested in a small amount of quartz wool to ensure retention of sample materials within the cup. The cups and quartz wool were ashed at 480 °C for 8 hours. Sample materials were loaded into cups within a Class 100 laminar flow isolation cabinet with HEPA filtered air to minimize the potential of particulate contamination.

During isotopic analysis, sample cups were evacuated and purged with helium prior to introduction into the EA. Reaction conditions were as follows: oxidation and reduction reactor temperatures were 1100 °C and 650 °C, respectively; helium carrier gas flow was 158 ml/min and the O₂ pulse was set for 45 seconds. Carbon dioxide generated during sample combustion was trapped within the EA in a molecular sieve trap. Following passage of the N₂ peak, the primary EA trap was heated and carbon dioxide was released to a secondary, silica gel-filled cryotrap which was immersed in liquid nitrogen. Trapping duration was calibrated using the EA thermal conductivity detector data to ensure complete collection of the CO₂ peak. Following collection of CO₂, the cryotrap gas flow was switched to a lower-flow He carrier gas (~1 mL/min) via an automated Vici Valco 6-port valve. The trap was warmed, and sample gas was released to the IRMS through an Agilent CarboBond capillary chromatography column (25 m x 0.53 mm x 5 µm). The resulting raw carbon isotope data are blank-corrected using direct blank subtraction and normalized to the VPDB scale using the two-point correction scheme⁷⁵ with the international reference materials NIST 1547-Peach Leaves (δ¹³C = 26.0 ± 0.2‰) and IAEA C6-sucrose (δ¹³C = −10.45 ± 0.03‰) which are run as solids. The reference materials USGS 61- Caffeine (−35.05‰) and USGS 62-Caffeine (−14.79‰) are dissolved in UV-treated MilliQ water, dispensed in known quantities, and have long-term laboratory reproducibility of ±0.47‰ for USGS 62 and ±0.39‰ for USGS 61 over a 20–90 nanomole range. Reproducibility of the δ¹³C values of reference materials with carbon contents greater than 50 nanomoles is ±0.3‰ (1 sd) and is equivalent to that previously reported⁷⁴.

TOF-SIMS analyses

The morphology and chemical composition of the fragments were characterized non-destructively by SEM and TOF-SIMS to confirm the presence of organic tissue well suited for carbon isotope analysis. Time-Of-Flight Secondary Ion Mass Spectrometry (TOF-SIMS) was used to produce semi-quantitative maps of the elemental composition in sample #15153 (Supplementary Fig. 5).

This technique bombards the surface with Bi ions that causes a collision cascade in the uppermost atom layers of the specimen (~10 nm). This releases secondary ions that are accelerated in an electric field and their time of flight to the detector in vacuum is a function of their mass and sample depth^76,77,78. Supplementary Fig. 5 shows elements bound to organic matter (incl. C, N, P) in the sample that produce polyatomic charged species, such as CN- and CNO- when emitted from the same sample depth. The correlations of organic-bound elements enable us to distinguish the presence of organic carbon from inorganic phases (e.g. carbonate minerals) in the sample.

Palaeoclimate modelling

The relatively fast coupled Earth-system model of intermediate complexity CLIMBER-3α was used to simulate the Devonian climate with an atmospheric CO₂ level of 500 ppm. CLIMBER-3α encompasses a modified version of the ocean circulation model (MOM3^79,80) with a horizontal resolution of 3.75° x 3.75° and 24 vertical levels, a dynamic/thermodynamic sea-ice model⁸¹ with the same resolution and a fast atmospheric model⁸² of 22.5° longitudinal and 7.5° latitudinal resolution. The model does not explicitly model ice sheet growth on the continents, but snow cover on the continents is considered. The model was run for Lower Devonian boundary conditions (415 Ma) in terms of continental configuration, solar luminosity and vegetation cover⁸³. Based on previous results⁶⁷, three different orbital configurations were explored: the standard configuration (obliquity 23.5°, eccentricity e = 0) as well as cold (obliquity 22.0°, eccentricity e = 0) and warm (obliquity 24.5°, eccentricity e = 0.069, precession angle 0°) orbital configurations (supplementary Fig. 17). A sensitivity analysis considering seasonal surface air temperatures and sea-ice distribution for these different insolation patterns is shown in Supplementary table 6. For determining climate variables at specific proxy locations (see Supplementary Data 2) the simulated values were bilinearly interpolated on a 1°-by-1° grid and evaluated using coordinates transformed from the present-day values using GPlates.

Long-term global biogeochemical modelling

We used the Carbon-Oxygen-Phosphorous-Sulphur Evolution (COPSE) model to predict the histories of atmospheric pCO₂, pO₂ and ocean composition over the Phanerozoic (550 Ma–today). This forward modelling approach enables hypothesis testing of mechanistic cause-effect relationships in the Earth system. A set of coupled differential equations describing the dynamic evolution of the C, O, P and S cycles were solved using an inbuilt variable timestep solver for ‘stiff’ Ordinary Differential Equation systems in Matlab®. The most recent version of the model (COPSE Reloaded, denoted ‘CR’; Supplementary Figs. 17–18) were adapted and modified^48,49. As input, we let the C/P ratio of buried terrestrial biomass increase with the colonization of land by non-vascular plants as in CR and adjusted the forcings on plant weathering (W), plant evolution (E), selective P weathering (F), volcanic outgassing (D) in a manner to simultaneously fit the effect of shallow vascular ecosystems on weathering processes and produce outputs consistent with the palaeorecords. Figure 3 shows the atmospheric pO₂ and pCO₂ trajectories predicted by the revised version of CR with forcings shown in supplementary Fig. 21. Further details on the revised COPSE modelling are found in the supplementary text section S5.

Data availability

All data are available in the main text or the supplementary materials.

Code availability

The models developed in this analysis are made available at University of Copenhagen’s Electronic Research Data Archive (ERD). https://doi.org/10.17894/ucph.214a6434-b7eb-4e62-aaac-afefc1247da4. The source code for the CLIMBER paleoclimate model used in this study is archived at the Potsdam Institute for Climate Impact Research and is made available upon request. The postprocessing scripts used to analyse the CLIMBER data and to generate the paleoclimate model figures presented in the study are stored on ERDA.

References

1. Stein, W. E. et al. Mid-devonian archaeopteris roots signal revolutionary change in earliest fossil forests. Curr. Biol. 30, 421–431.e2 (2020).Article CAS Google Scholar 
2. Foster, G. L., Royer, D. L. & Lunt, D. J. Future climate forcing potentially without precedent in the last 420 million years. Nat. Commun. 8, 14845 (2017).Article ADS CAS Google Scholar 
3. Berner, R. A. Paleozoic atmospheric CO₂: importance of solar radiation and plant evolution. Science 261, 68–70 (1993).Article ADS CAS Google Scholar 
4. Walker, J. C. G., Hays, P. B. & Kasting, J. F. A negative feedback mechanism for the long-term stabilization of Earth’s surface temperature. J. Geophys. Res. 86, 9776 (1981).Article ADS CAS Google Scholar 
5. D’Antonio, M., Ibarra, D. E. & Boyce, C. K. Land plant evolution decreased, rather than increased, weathering rates. Geology https://doi.org/10.1130/G46776.1 (2019).Article Google Scholar 
6. Caves, J. K., Jost, A. B., Lau, K. V. & Maher, K. Cenozoic carbon cycle imbalances and a variable weathering feedback. Earth Planet. Sci. Lett. 450, 152–163 (2016).Article ADS CAS Google Scholar 
7. Maher, K. & Chamberlain, C. P. Hydrologic regulation of chemical weathering and the geologic carbon cycle. Science 343, 1502–1504 (2014).Article ADS CAS Google Scholar 
8. Goddéris, Y. et al. Onset and ending of the late Palaeozoic ice age triggered by tectonically paced rock weathering. Nat. Geosci. 10, 382–386 (2017).Article ADS Google Scholar 
9. Ibarra, D. E. et al. Modeling the consequences of land plant evolution on silicate weathering. Am. J. Sci. 319, 1–43 (2019).Article ADS CAS Google Scholar 
10. Berner, R. A. Atmospheric carbon dioxide levels over phanerozoic time. Science 249, 1382–1386 (1990).Article ADS CAS Google Scholar 
11. Berner, R. A. GEOCARBSULF. A combined model for Phanerozoic atmospheric O₂ and CO₂. Geochimica et Cosmochimica Acta 70, 5653–5664, (2006).
12. Lenton, T. M. et al. Earliest land plants created modern levels of atmospheric oxygen. Proc. Natl. Acad. Sci. USA 113, 9704–9709 (2016).Article ADS CAS Google Scholar 
13. Ghienne, J.-F. et al. A Cenozoic-style scenario for the end-Ordovician glaciation. Nat. Commun. 5, 4485 (2014).Article ADS CAS Google Scholar 
14. Pohl, A. et al. Glacial onset predated Late Ordovician climate cooling. Paleoceanography 31, 800–821 (2016).Article ADS Google Scholar 
15. Royer, D. L. Stomatal density and stomatal index as indicators of paleoatmospheric CO₂ concentration. Rev. Palaeobot. Palynol. 114, 1–28 (2001).Article Google Scholar 
16. Mora, C. I., Driese, S. G. & Colarusso, L. A.Middleto Late Paleozoic Atmospheric CO$_2$ Levels from Soil Carbonate and Organic Matter. Sci. New Ser. 271, 1105–1107 (1996).CAS Google Scholar 
17. Ekart, D. D. A 400 million year carbon isotope record of pedogenic carbonate; implications for paleoatomospheric carbon dioxide. Am. J. Sci. 299, 805–827 (1999).Article ADS CAS Google Scholar 
18. Witkowski, C. R., Weijers, J. W. H., Blais, B., Schouten, S. & Sinninghe Damsté, J. S. Molecular fossils from phytoplankton reveal secular P co ₂ trend over the Phanerozoic. Sci. Adv. 4, eaat4556 (2018).Article ADS CAS Google Scholar 
19. Yapp, C. J. & Poths, H. Ancient atmospheric CO₂ pressures inferred from natural goethites. Nature 355, 342–344 (1992).Article ADS CAS Google Scholar 
20. McElwain, J. & Chaloner, W. Stomatal density and index of fossil plants track atmospheric carbon dioxide in the palaeozoic. Ann. Bot. 76, 389–395 (1995).Article Google Scholar 
21. Yapp, C. J. A possible goethite-iron(III) carbonate solid solution and the determination of CO₂ partial pressures in low-temperature geologic systems. Chem. Geol. 64, 259–268 (1987).Article ADS CAS Google Scholar 
22. Russell, J. D., Paterson, E., Fraser, A. R. & Farmer, V. C. Adsorption of carbon dioxide on goethite (α-FeOOH) surfaces, and its implications for anion adsorption. J. Chem. Soc. Faraday Trans. 1: Phys. Chem. Condens. Phases 71, 1623 (1975).Article CAS Google Scholar 
23. Breecker, D. O., Sharp, Z. D. & McFadden, L. D. Atmospheric CO₂ concentrations during ancient greenhouse climates were similar to those predicted for A.D. 2100. Proc. Natl. Acad. Sci. 107, 576–580 (2010).Article ADS CAS Google Scholar 
24. Edwards, D. Stomata in early land plants: an anatomical and ecophysiological approach. J. Exp. Bot. 49, 255–278 (1998).Article Google Scholar 
25. Witkowski, C. R., van der Meer, M. T. J., Smit, N. T., Sinninghe Damsté, J. S. & Schouten, S. Testing algal-based pCO₂ proxies at a modern CO₂ seep (Vulcano, Italy). Sci. Rep. 10, 10508 (2020).Article ADS CAS Google Scholar 
26. Freeman, K. H. & Hayes, J. M. Fractionation of carbon isotopes by phytoplankton and estimates of ancient CO ₂ levels. Glob. Biogeochem. Cycles 6, 185–198 (1992).Article ADS CAS Google Scholar 
27. Franks, P. J. et al. New constraints on atmospheric CO₂ concentration for the Phanerozoic. Geophys. Res. Lett. 41, 4685–4694 (2014).Article ADS CAS Google Scholar 
28. Diefendorf, A. F., Mueller, K. E., Wing, S. L., Koch, P. L. & Freeman, K. H. Global patterns in leaf 13C discrimination and implications for studies of past and future climate. Proc. Natl Acad. Sci. 107, 5738–5743 (2010).Article ADS CAS Google Scholar 
29. Schubert, B. A. & Jahren, A. H. Global increase in plant carbon isotope fractionation following the Last Glacial Maximum caused by increase in atmospheric pCO₂. Geology 43, 435–438 (2015).Article ADS CAS Google Scholar 
30. Schubert, B. A. & Jahren, A. H. The effect of atmospheric CO₂ concentration on carbon isotope fractionation in C3 land plants. Geochimica et. Cosmochimica Acta 96, 29–43 (2012).Article ADS CAS Google Scholar 
31. Cui, Y. & Schubert, B. A. A 23 m.y. record of low atmospheric CO₂. Geology 48, 888–892 (2020).Article ADS CAS Google Scholar 
32. Diefendorf, A. F., Leslie, A. B. & Wing, S. L. Leaf wax composition and carbon isotopes vary among major conifer groups. Geochimica et. Cosmochimica Acta 170, 145–156 (2015).Article ADS CAS Google Scholar 
33. Lomax, B. H., Lake, J. A., Leng, M. J. & Jardine, P. E. An experimental evaluation of the use of Δ13C as a proxy for palaeoatmospheric CO₂. Geochimica et. Cosmochimica Acta 247, 162–174 (2019).Article ADS CAS Google Scholar 
34. Stein, R. A., Sheldon, N. D. & Smith, S. Rapid response to anthropogenic climate change by Thuja occidentalis: implications for past climate reconstructions and future climate predictions. PeerJ 7, e7378 (2019).Article Google Scholar 
35. Stein, R. A., Sheldon, N. D. & Smith, S. Y. C3 plant carbon isotope discrimination does not respond to CO₂ concentration on decadal to centennial timescales. N. Phytologist 229, 2576–2585 (2021).Article CAS Google Scholar 
36. Jardine, P. E. & Lomax, B. H. A 23 m.y. record of low atmospheric CO₂. Geology 49, e523 (2021).Article ADS Google Scholar 
37. Franks, P. J. & Farquhar, G. D. The mechanical diversity of stomata and its significance in gas-exchange control. Plant Physiol. 143, 78–87 (2007).Article CAS Google Scholar 
38. Kohn, M. J. Carbon isotope compositions of terrestrial C3 plants as indicators of (paleo)ecology and (paleo)climate. Proc. Natl Acad. Sci. 107, 19691–19695 (2010).Article ADS CAS Google Scholar 
39. Fletcher, B. J., Brentnall, S. J., Quick, W. P. & Beerling, D. J. BRYOCARB: A process-based model of thallose liverwort carbon isotope fractionation in response to CO₂, O₂, light and temperature. Geochimica et. Cosmochimica Acta 70, 5676–5691 (2006).Article ADS CAS Google Scholar 
40. Porter, A. S., Yiotis, C., Montañez, I. P. & McElwain, J. C. Evolutionary differences in Δ13C detected between spore and seed bearing plants following exposure to a range of atmospheric O₂:CO₂ ratios; implications for paleoatmosphere reconstruction. Geochimica et. Cosmochimica Acta 213, 517–533 (2017).Article ADS CAS Google Scholar 
41. Porter, A. S. et al. Testing the accuracy of new paleoatmospheric CO₂ proxies based on plant stable carbon isotopic composition and stomatal traits in a range of simulated paleoatmospheric O₂:CO₂ ratios. Geochimica et. Cosmochimica Acta 2019, 69–90 (2019).Article ADS Google Scholar 
42. Boyce, C. K. How green was Cooksonia? The importance of size in understanding the early evolution of physiology in the vascular plant lineage. Paleobiology 34, 179–194 (2008).Article Google Scholar 
43. Streel, M. Late Frasnian–Famennian climates based on palynomorph analyses and the question of the Late Devonian glaciations. Earth-Sci. Rev. 52, 121–173 (2000).Article ADS Google Scholar 
44. Cleal, C. & Cascales-Miñana, B. Evolutionary Floras – revealing large-scale patterns in Palaeozoic vegetation history. J. Palaeosci. 28, 1–20 (2021).
45. Dahl, T. W. & Arens, S. K. M. The impacts of land plant evolution on Earth’s climate and oxygenation state – An interdisciplinary review. Chem. Geol. 547, 119665 (2020).Article ADS CAS Google Scholar 
46. McMahon, W. J. & Davies, N. S. Evolution of alluvial mudrock forced by early land plants. Science 359, 1022–1024 (2018).Article ADS CAS Google Scholar 
47. Libertín, M., Kvaček, J., Bek, J., Žárský, V. & Štorch, P. Sporophytes of polysporangiate land plants from the early Silurian period may have been photosynthetically autonomous. Nat. Plants 4, 269–271 (2018).Article Google Scholar 
48. Lenton, T. M., Daines, S. J. & Mills, B. J. W. COPSE reloaded: An improved model of biogeochemical cycling over Phanerozoic time. Earth-Sci. Rev. 178, 1–28 (2018).Article ADS CAS Google Scholar 
49. Krause, A. J. et al. Stepwise oxygenation of the Paleozoic atmosphere. Nat. Commun. 9, 4081 (2018).Article ADS Google Scholar 
50. Porada, P. et al. High potential for weathering and climate effects of non-vascular vegetation in the Late Ordovician. Nat. Commun. 7, 12113 (2016).Article ADS CAS Google Scholar 
51. Halder, S., Arens, S. K. M., Jensen, K., Dahl, T. W. & Porada, P. A dynamic local-scale vegetation model for lycopsids (LYCOm v1.0). Geosci. Model Dev. 15, 2325–2343 (2022).Article ADS CAS Google Scholar 
52. Bergman, N. M. COPSE: A new model of biogeochemical cycling over Phanerozoic time. Am. J. Sci. 304, 397–437 (2004).Article ADS CAS Google Scholar 
53. McKenzie, N. R. et al. Continental arc volcanism as the principal driver of icehouse-greenhouse variability. Science 352, 444–447 (2016).Article ADS CAS Google Scholar 
54. Glasspool, I. J., Edwards, D. & Axe, L. Charcoal in the Silurian as evidence for the earliest wildfire. Geol. 32, 381 (2004).Article ADS Google Scholar 
55. Glasspool, I. J. & Scott, A. C. Phanerozoic concentrations of atmospheric oxygen reconstructed from sedimentary charcoal. Nat. Geosci. 3, 627–630 (2010).Article ADS CAS Google Scholar 
56. Sønderholm, F. & Bjerrum, C. J. Minimum levels of atmospheric oxygen from fossil tree roots imply new plant−oxygen feedback. Geobiology 19, 250–260 (2021).Article Google Scholar 
57. Algeo, T. J., Scheckler, S. E. & Scott, A. C. Terrestrial-marine teleconnections in the devonian: links between the evolution of land plants, weathering processes, and marine anoxic events [and discussion]. Philos. Trans.: Biol. Sci. 353, 113–130 (1998).Article Google Scholar 
58. Montañez, I. P. et al. Climate, pCO₂ and terrestrial carbon cycle linkages during late Palaeozoic glacial–interglacial cycles. Nat. Geosci. 9, 824–828 (2016).Article ADS Google Scholar 
59. Arora, V. K. & Montenegro, A. Small temperature benefits provided by realistic afforestation efforts. Nat. Geosci. 4, 514–518 (2011).Article ADS CAS Google Scholar 
60. Grossman, E. L. & Joachimski, M. M. Oxygen Isotope Stratigraphy. in Geologic Time Scale 2020 279–307 (Elsevier, 2020). https://doi.org/10.1016/B978-0-12-824360-2.00010-3.
61. Chen, B. et al. Devonian paleoclimate and its drivers: A reassessment based on a new conodont δ18O record from South China. Earth-Sci. Rev. 222, 103814 (2021).Article CAS Google Scholar 
62. van Geldern, R. et al. Carbon, oxygen and strontium isotope records of Devonian brachiopod shell calcite. Palaeogeogr., Palaeoclimatol., Palaeoecol. 240, 47–67 (2006).Article Google Scholar 
63. Liu, Y., Peltier, W. R., Yang, J. & Vettoretti, G. The initiation of Neoproterozoic “snowball” climates in CCSM3: the influence of paleocontinental configuration. Clim. Past 23, 2555–2577 (2013).
64. Feulner, G. Climate simulations of Neoproterozoic snowball Earth events: Similar critical carbon dioxide levels for the Sturtian and Marinoan glaciations. Earth Planet. Sci. Lett. 404, 200–205 (2014).Article ADS CAS Google Scholar 
65. Feulner, G. Formation of most of our coal brought Earth close to global glaciation. Proc. Natl. Acad. Sci. USA 114, 11333–11337 (2017).Article ADS CAS Google Scholar 
66. De Vleeschouwer, D., Crucifix, M., Bounceur, N. & Claeys, P. The impact of astronomical forcing on the Late Devonian greenhouse climate. Glob. Planet. Change 120, 65–80 (2014).Article ADS Google Scholar 
67. Brugger, J., Hofmann, M., Petri, S. & Feulner, G. On the sensitivity of the devonian climate to continental configuration, vegetation cover, orbital configuration, CO₂ concentration, and insolation. Paleoceanogr. Paleoclimatology 34, 1375–1398 (2019).Article ADS Google Scholar 
68. Isbell, J. L. et al. Glacial paradoxes during the late Paleozoic ice age: Evaluating the equilibrium line altitude as a control on glaciation. Gondwana Res. 22, 1–19 (2012).Article ADS Google Scholar 
69. Gindre, L., Le Heron, D. & Bjørnseth, H. M. High resolution facies analysis and sequence stratigraphy of the Siluro-Devonian succession of Al Kufrah basin (SE Libya). J. Afr. Earth Sci. 76, 8–26 (2012).Article ADS Google Scholar 
70. Gerrienne, P., Bergamaschi, S., Pereira, E., Rodrigues, M.-A. C. & Steemans, P. An Early Devonian flora, including Cooksonia, from the Paraná Basin (Brazil). Rev. Palaeobot. Palynol. 116, 19–38 (2001).Article Google Scholar 
71. Cookson, I. IV—On plant-remains from the Silurian of Victoria, Australia, that extend and connect floras hitherto described. Philos. Trans. R. Soc. Lond. B 225, 127–148 (1935).Article ADS Google Scholar 
72. Earp, C. Costulatotheca schleigeri (Hyolitha: Orthothecida) from the Walhalla Group (Early Devonian) at Mount Pleasant, central Victoria, Australia. Alcheringa: Australas. J. Palaeontol. 43, 220–227 (2019).Article Google Scholar 
73. Cohen, K. M., Finney, S. C., Gibbard, P. L. & Fan, J.-X. The ICS international chronostratigraphic chart. Episodes 36, 199–204 (2013).Article Google Scholar 
74. Polissar, P. J., Fulton, J. M., Junium, C. K., Turich, C. C. & Freeman, K. H. Measurement of ¹³ C and ¹⁵ N Isotopic Composition on Nanomolar Quantities of C and N. Anal. Chem. 81, 755–763 (2009).Article CAS Google Scholar 
75. Coplen, T. B. et al. New guidelines for δ ¹³ C measurements. Anal. Chem. 78, 2439–2441 (2006).Article CAS Google Scholar 
76. Bertrand, P. & Lu-Tao, W. Time-of-Flight Secondary Ion Mass Spectrometry (ToF-SIMS). in Microbeam and Nanobeam Analysis (eds. Benoit, D., Bresse, J.-F., Van’t Dack, L., Werner, H. & Wernisch, J.) 167–182 (Springer Vienna, 1996). https://doi.org/10.1007/978-3-7091-6555-3_8.
77. Belu, A. M., Graham, D. J. & Castner, D. G. Time-of-flight secondary ion mass spectrometry: techniques and applications for the characterization of biomaterial surfaces. Biomaterials 24, 3635–3653 (2003).Article CAS Google Scholar 
78. Sodhi, R. N. S. Time-of-flight secondary ion mass spectrometry (TOF-SIMS):—versatility in chemical and imaging surface analysis. Analyst 129, 483–487 (2004).Article ADS CAS Google Scholar 
79. Pacanowski, R. & Griffies, S. The MOM 3 manual, alpha version. NOAA. Geophysical Fluid Dynamics Laboratory (1999).
80. Hofmann, M. & M.A., Morales Maqueda. Performance of a second-order moments advection scheme in an Ocean General Circulation Model. J. Geophys. Res. 111, C05006 (2006).ADS Google Scholar 
81. Fichefet, T. & Morales Maqueda, M. A. Sensitivity of a global sea ice model to the treatment of ice thermodynamics and dynamics. J. Geophys. Res. 102, 12609–12646 (1997).Article ADS Google Scholar 
82. Petoukhov, V. et al. CLIMBER-2: a climate system model of intermediate complexity. Part I: model description and performance for present climate. Clim. Dyn. 16, 1–17 (2000).Article Google Scholar 
83. Scotese, C. R. Atlas of Devonian Paleogeographic Maps, PALEOMAP Atlas for ArcGIS, The Late Paleozoic, Maps 65-72, Mollweide Projection. (2014).
84. Belcher, C. M. & McElwain, J. C. Limits for combustion in low O₂ redefine paleoatmospheric predictions for the mesozoic. Science 321, 1197–1200 (2008).Article ADS CAS Google Scholar 

Download references

Acknowledgements

We thank T. Ziegler and F. McSweeney for providing specimens from the palaeobotanical collection of Museum Victoria, O. Seberg and A. A. Pedersen for providing access to lycophyte cultures in the greenhouse of the Botanical Garden, University of Copenhagen, and to F. Sønderholm and C. J. Bjerrum for assistance and use of pCO₂ metre. T.W.D. was funded by the Carlsberg Foundation through its Distinguished Associate Professsor program (grant no. CF16–0876) and the Danish Council for Independent Research (grant nos. 7014-00295B, 8102-00005B). C. J. was funded by the National Science Foundation (NSF EAR-1455258). The authors gratefully acknowledge the European Regional Development Fund (ERDF), the German Federal Ministry of Education and Research and the Land Brandenburg for supporting this project by providing resources on the high-performance computer system at the Potsdam Institute for Climate Impact Research. J. B. was funded through the VeWA consortium (Past Warm Periods as Natural Analogues of our high-CO₂ Climate Future) by the LOEWE programme of the Hessen Ministry of Higher Education, Research and the Arts, Germany. B.H.L acknowledges funding from the NERC (grant nos NE/R001324/1, NE/T00392/1).

Author information

Authors and Affiliations

1. Globe institute, University of Copenhagen; Øster Voldgade 5−7, Copenhagen, DenmarkTais W. Dahl & Magnus A. R. Harding
2. State Key Laboratory of Geological Processes and Mineral Resources, China University of Geosciences, Wuhan, ChinaTais W. Dahl
3. Sino-Danish College (SDC), University of Chinese Academy of Sciences, Beijing, ChinaMagnus A. R. Harding
4. Earth System Analysis, Potsdam Institute for Climate Impact Research, Member of the Leibniz Association, Potsdam, GermanyJulia Brugger & Georg Feulner
5. Senckenberg Biodiversity and Climate Research Centre, Frankfurt am Main, GermanyJulia Brugger
6. Center for Integrative Petroleum Research, King Fahd University of Petroleum and Minerals, Dhahran, Saudi ArabiaKion Norrman
7. School of Biosciences, University of Nottingham, Sutton Bonington Campus, Leicestershire, UKBarry H. Lomax
8. Department of Earth and Environmental Sciences, Syracuse University; Syracuse, New York, USAChristopher K. Junium

Contributions

TWD designed the research with BHL’s assistance; C.K.J., M.A.R.H., K.N., and T.W.D collected data; G.F., J.B., B.H.L., and T.W.D. performed research; T.W.D., C.K.J., and G.F. acquired funding. T.W.D. wrote the paper with input from all authors.

Corresponding author

Correspondence to Tais W. Dahl.

Ethics declarations

Competing interests

The authors declare no competing interests

Peer review

Peer review information

Nature Communications thanks Alexander Krause and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Supplementary information

Peer Review File

Description of Additional Supplementary Files

Supplementary Dataset 1

Supplementary Dataset 2

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and Permissions

About this article

Cite this article

Dahl, T.W., Harding, M.A.R., Brugger, J. et al. Low atmospheric CO₂ levels before the rise of forested ecosystems. Nat Commun 13, 7616 (2022). https://doi.org/10.1038/s41467-022-35085-9

Download citation

• Received18 January 2022
• Accepted18 November 2022
• Published20 December 2022
• DOIhttps://doi.org/10.1038/s41467-022-35085-9