Permo–Triassic boundary carbon and mercury cycling linked to terrestrial ecosystem collapse

Records suggest that the Permo–Triassic mass extinction (PTME) involved one of the most severe terrestrial ecosystem collapses of the Phanerozoic. However, it has proved difficult to constrain the extent of the primary productivity loss on land, hindering our understanding of the effects on global biogeochemistry. We build a new biogeochemical model that couples the global Hg and C cycles to evaluate the distinct terrestrial contribution to atmosphere–ocean biogeochemistry separated from coeval volcanic fluxes. We show that the large short-lived Hg spike, and nadirs in δ202Hg and δ13C values at the marine PTME are best explained by a sudden, massive pulse of terrestrial biomass oxidation, while volcanism remains an adequate explanation for the longer-term geochemical changes. Our modelling shows that a massive collapse of terrestrial ecosystems linked to volcanism-driven environmental change triggered significant biogeochemical changes, and cascaded organic matter, nutrients, Hg and other organically-bound species into the marine system. The environmental changes at the Permian–Triassic boundary are thought to have been caused primarily by volcanic eruptions. Here the authors develop a model to show that the loss of ecosystems on land and consequent massive terrestrial biomass oxidation triggered large biogeochemical changes in the oceans at the time of the marine mass extinction.

T he Permo-Triassic mass extinction (PTME) is the largest known extinction in Earth′s history, with the loss of~90% of species in the sea and~70% of species on land [1][2][3][4] . The PTME has been causally linked to the emplacement of the Siberian Traps Large Igneous Province (LIP) and associated volcanic gas emissions (especially CO 2 , SO 2 and halogens), via widespread environmental changes such as warming and oceanic anoxia [5][6][7][8] . The PTME also saw a crisis in terrestrial ecosystems, with loss of plant diversity, increased wildfire activity and consequent enhanced soil erosion [9][10][11][12][13][14] . Recent work has shown that the disruption of vegetation started before and culminated at the marine extinction level (Fig. 1), implying that the environmental disaster impacted terrestrial ecosystems first 12,14 . The cause of the terrestrial mass extinction is still unclear, and several kill mechanisms have been hypothesised. For example, a shift from a humid warm climate to an unstable highly seasonal climate and an associated increase in wildfires affected the equatorial Permo-Triassic peatlands, drastically reducing the abundance and diversity of the flora 14 ; abnormal pollen and spores found in different localities around the world during the PTME interval suggest widespread mutagenesis possibly linked to an increase in UV-B radiation due to ozone depletion 15,16 ; a terrestrial S-isotope record from the Karoo basin in South Africa could indicate volcanically driven acid rain at the P-T transition 17 that might have also severely impacted the flora. Whilst the taxonomic losses in terrestrial ecosystems are becoming clearer 12,18 , and local, enhanced input of terrestrial material into marine environments has been recorded 9,11 , the biogeochemical impacts and feedbacks on the exogenic C cycle are not known. It is possible that these impacts were severe: the PTME represents the largest, and maybe the only, known mass extinction of insects 19 , suggesting that there may have been a substantial decrease in available food sources at the lowest levels of the food chain.
The PTME is marked by an approximately two-to four-fold increase in marine sedimentary Hg concentration with respect to background levels during a~400 kyr interval also characterised by negative δ 13 C values 20 , which implies a relatively long-term injection of Hg and 13 C-depleted CO 2 into the atmosphereland-ocean system during this time. Superimposed on this trend is a prominent, short-lived Hg spike, which is usually expressed as Hg/TOC, given the affinity of Hg with organic matter, and which is coincident with the collapse of the terrestrial ecosystems 14 , the onset of the marine mass extinction interval and a sharp minimum in δ 13 C values (Fig. 1). While higher Hg/TOC values have been reported preceding the marine extinction in the deep-water settings of Japan 21 , these are an artefact of normalisation to values of TOC which are below the analytical detection limit (<<0.1%), and do not track increases in Hg concentrations (see TOC and Hg data in the Supplementary Information of ref. 21 ). The Hg record is interpreted as evidence of increased Hg input into the Earth's surface system from the Siberian Traps 13,21-24 . However, Hg is also stored in large quantities in terrestrial biomass and soils [25][26][27] . The mobilisation of these terrestrial reservoirs during the PTME could also result in the increased loading of Hg to aquatic environments, even without an elevated volcanic Hg flux to the atmosphere 23,[27][28][29] . Any such changes in the soil and biomass carbon reservoir would also have direct implications for the release of C to the atmosphere and the sedimentary δ 13 C record.
Hg isotopes can be used to better understand how Hg has been transported into the sedimentary environment. Mass-independent fractionation (MIF-denoted by Δ) occurs due to aqueous photoreduction of Hg 2+ to Hg 0 ðgÞ that takes place in surface waters and  Fig. 1 Permian-Triassic events. Existing records show the terrestrial extinction started before, and its collapse culminated at, the marine extinction interval, where large Hg and Hg/TOC spikes, negative δ 202 Hg, negative Δ 199 Hg and more negative δ 13 C values (EP. I 30 ) are recorded. After this event, the δ 13 C record shows a second minimum (Ep. II), while Hg concentration decreases, but remains relatively higher than background values, and δ 202 Hg and Δ 199 Hg rebound towards more positive values. This pattern is seen in many P-T successions, in terrestrial and shallow-to deep-water settings. Meishan section: U-Pb ages from Burgess et al. 34 ; TOC and C isotopes (δ 13 C) from Xie et al. 30 ; Hg concentrations, δ 202 Hg and Δ 199 Hg data are from Grasby et al. 23 and Wang et al. 20 . Chinahe section: all data and correlation to Meishan section are from Chu et al. 14 . Note that to plot Hg (ppb) and Hg/TOC data, the same scale is used. MDF mass-dependent fractionation, MIF mass-independent fractionation, EP. episode. in clouds 29 . This results in positive Δ 199 Hg values in the remaining water Hg 2+ pool, hence positive Δ 199 Hg values in sediments dominated by atmospheric Hg 2+ deposition. Conversely, slightly negative Δ 199 Hg values characterise the terrestrial reservoir (soil and biomass), which primarily captures Hg 0 ðgÞ 29 . During Hg uptake by plants additional mass-dependent fractionation (MDF) and MIF occur, resulting in more negative δ 202 Hg and Δ 199 Hg values 29 . The Permian-Triassic boundary shallow-water record of Meishan shows a prominent negative δ 202 Hg excursion in correspondence to the Hg spike coupled to a small negative Δ 199 Hg shift, but deeper water successions show persistent positive Δ 199 Hg 20,21,23,24 . Therefore, isotope data appear to indicate that the Hg was transported to the deep-water settings mostly via the atmosphere, and to the shallow-water settings via continental runoff and via the atmosphere 20,21,23,24 .
The δ 13 C records at the Permian-Triassic boundary show two minima 30,31 , which are here called EP. I and II (EP. = episode) following ref. 30 (Fig. 1). The observed δ 13 C trends are similarly recorded in different depositional settings and by different substrates (carbonates, bulk organic matter, separate plant remains) 14,[30][31][32] , strongly indicating that they represent actual changes in the C-isotope composition of the reservoirs of the exogenic C cycle. The Hg concentration spike occurs in the same interval as the minimum in δ 13 C values associated with the PTME (EP. I, Fig. 1). At Meishan, where the chronostratigraphic framework is well established 33,34 , the initial large Hg spike occurs about 60 Kyrs after the onset of the carbonisotope perturbation (Fig. 1). Data from non-marine end-Permian successions confirms this diachrony 13,14 (Fig. 1) and show that the Hg spike is also coincident with a sudden decrease in Total Organic Carbon (TOC) values to almost zero 14 (Fig. 1). The combination of geochemical and palaeontological data from these sections shows that the terrestrial ecosystem disruption started with the onset of the carbon-isotope perturbation and climaxed at the very sharp δ 13 C minimum (EP. I), coincident with the Hg spike ( Fig. 1), and the start of the main marine extinction interval. This mismatch in both timings and fluxes between the C and Hg cycles at the PTME, suggest that the 13 C-depleted C and the Hg came from multiple sources.
Overall, the interplay between volcanism and terrestrial reservoir changes in controlling PTME biogeochemistry is not well known, and previous attempts to model the δ 13 C record have been fundamentally hampered because of the lack of an independent tracer of the C source. To overcome this, we use published records of δ 13 C and Hg systematics to jointly constrain a new coupled C and Hg biogeochemical model. Our model shows that the large, sudden geochemical shifts at the PTME are best explained by a massive pulse of terrestrial biomass oxidation, while Siberian Traps volcanism can explain the longer-term geochemical changes.

Results and discussion
A coupled Hg-C cycle model. Figure 2 shows the biogeochemical box model developed here. The full model derivation follows in the 'Methods' section. The model combines a multi-box sediment-ocean-atmosphere carbon-alkalinity cycle (based on previous work [35][36][37], with the global mercury cycle 38,39 . The ocean is split into 'surface', 'high-latitude' and 'deep' boxes. It considers ocean circulation and carbonate speciation, and contains a simplified organic carbon cycle in which burial rates are prescribed. As well as computing the global C and Hg cycles it also computes δ 13 C of all C reservoirs and δ 202 Hg of the ocean reservoirs. A full atmosphere-ocean model of δ 202 Hg would require dynamic biosphere reservoirs, which would greatly increase model complexity.
We therefore simplify the system to a mixing model for marine δ 202 Hg, in which atmospheric and riverine inputs have different isotopic signatures. The atmosphere is assumed to have δ 202 Hg atm ¼ À1 ‰, and riverine input is assumed to have δ 202 Hg runoff ¼ À3 ‰, following ref. 40 .
The model is set up for the late Permian by reducing the solar constant to that of 250 Ma, and increasing the background tectonic CO 2 degassing rate to 1.5 times the present day (D = 1.5), in line with estimates for the Late Permian 41 . To obtain the observed pre-event ocean-atmosphere δ 13 C composition of 3.5‰, we set the rate of land-derived organic carbon burial to 20% higher than present day, and adjust the composition of the weathered carbonate reservoir to 3‰. This is consistent with high terrestrial productivity and C burial in the Permian (e.g., coal forests and mires) and rapid recycling of more recently buried and 13 C-enriched carbonate material.
In the following paragraphs, we test two model end-member scenarios: (I) the release of volcanic and thermogenic Hg and C from Siberian Traps activity alone, and (II) with the additional release of Hg and C as a consequence of the collapse of the terrestrial ecosystems.
Volcanic and thermogenic degassing. We first model the release of volcanic/volcanogenic Hg and C from the Siberian Traps. Existing radioisotope data show that the extinction, the negative δ 13 C excursion and Hg spike might have all occurred during the intrusive phase of the Siberian Traps 20,34,42 . It is suggested that the emplacement of large sills caused the combustion or thermal decomposition of organic-rich sediments with the consequent release of thermogenic volatiles, such as C and Hg 20,22 . It has been proposed that over a~400 Kyrs intrusive phase the Siberian Traps emitted~7600-13,000 Mg yr −1 of volcanic Hg, which included both magmatic and coal-derived Hg 20,22,23 . Relating this Hg release to the background volcanic source is not straightforward because estimates of the background source vary, but taking the most likely present day range 25 (~90-360 Mg yr −1 ), and further constraining this by taking into account the need to balance overall sedimentary burial of Hg (~190 Mg yr −1 ), and the overall~50% increase in tectonic degassing in the late Palaeozoic relative to today 43 , we arrive at a best guess for the background late Permian Hg flux of 300 Mg yr −1 . This means that the Siberian Traps eruption increased the geogenic Hg input by a factor of~25-43 over~400 Kyrs.
To test this scenario, we model Siberian Traps intrusion by increasing the volcanic Hg source by 25-43-fold for 450 Kyrs, while also increasing the CO 2 source in line with estimates 44,45 for Siberian Traps degassing based on magma volumes and sediment intrusion (by 4-8 × 10 12 moles/year). The CO 2 released by contact metamorphism at the PTME is assumed to have an average δ 13 C composition of −25‰ 44 . Specifically, the input functions are: Here the first vector is time in millions of years and the second is the flux alteration at that time. Here CO 2ramp is the additional CO 2 release in mol yr −1 , and Hg ramp is the relative Hg degassing rate increase. For the duration of these pulses, the thermohaline circulation is also assumed to collapse due to warming and freshwater input 46 . We reduce the circulation term to 1 Sv over this period, which allows more rapid change in the model surface ocean C isotopes and Hg loading. This is a large reduction, and also reflects the simple structure of the model in which the entire low-latitude surface ocean is represented by a single box, and so is well-mixed. Figure 3a-d shows that this magnitude and timing of release of C and Hg is capable of driving the longer-term decline in carbonate δ 13 C, and the coeval long-term approximately twoto four-fold enrichment in shallow sediment Hg/TOC that is observed in Meishan. However, the model scenario does not capture the spike in Hg concentration, or nadir in δ 13 C (EP. I 30 in Fig. 1) that are coincident with the final stage of the terrestrial extinction. It also does not produce any substantial change in marine δ 202 Hg isotopes (Fig. 1d), because the primary Hg source to the ocean is the atmosphere for the full model run.
Within the model, we have also explored a scenario wherein the large Hg pulse represents a further rapid pulse of LIP volcanism. We have attempted this scenario in Supplementary Note 1 (Scenario I-2), where a 1 Kyr volcanic pulse is assumed to raise the Hg and C input rates by a further factor of 5. While the Hg/TOC can indeed be explained by an additional short-lived pulse of Hg, we require the total release rate of Hg to be~200 times greater than background levels, and even then, this scenario fails to reproduce any of the Hg isotope signature or the nadir in carbonate δ 13 C (Supplementary Fig. 1).
Terrestrial ecosystem collapse. For scenario II, we explore the additional effects of a geologically rapid (~1 Kyr) pulse of Hg and C as the result of the collapse of terrestrial ecosystems at the PTME. The magnitude of this Hg flux is again difficult to quantify precisely, and we explore an increase of 100-fold over background conditions. This level of increase represents the magnitude required to drive the sedimentary Hg signal that we observe, and is compatible with the available terrestrial biosphere Hg reservoir: total soil Hg is estimated to be on the order of~10 6 Mg Hg when considering a soil depth of~15 cm 47 . So, our model Hg delivery flux would require decimetre-scale soil organic matter oxidation over 1000 years, coincident with the PTME and the sharp EP. I negative δ 13 C shift 11 . The Hg pulse is delivered directly to the lowlatitude surface ocean via runoff in the model, and is accompanied by a pulse of 'soil oxidation' C which we assume raises the global rate of oxidative weathering by a factor of 30-a number chosen to have the observed level of impact on the C-isotope record while being compatible with the Hg input change. We also assume a cessation of terrestrial organic C burial. Terrestrial Hg deposition and erosion is not altered during the pulse as the fluxes are minor by comparison. The new model functions applied in addition to the longer-term inputs of scenario I are: Here, the first vector is time in millions of years, and the second shows flux multipliers at these times. C ramp ; O ramp and Hg bio denote the relative rate of land organic C burial, oxidative weathering and Hg runoff, respectively, and are set at 0, 30 and 100, respectively, for the duration of the 1-kyr pulse.
This 'biosphere' pulse causes a short-term large concentration spike in the shallow marine Hg reservoir and its sediments, which is superimposed on the volcanically driven changes (Fig. 3e-h). The Hg spike is far larger than would be expected from simply increasing the volcanic source by the same amount because the biospheric Hg is delivered directly to the surface ocean and sedimentation occurs mostly on the shelf. With the inclusion of terrestrial C oxidation and cessation of terrestrial carbon burial, the model also replicates the transient shift to more negative δ 13 C values recorded at the marine extinction interval (EP. I 30,31 in Fig. 1): Terrestrial C oxidation is a source of isotopically light C 48 . The model now also shows a sharp negative δ 202 Hg shift in the shallow ocean box, which is triggered by increased Hg riverine input, but shows no change in the deeper ocean, where the source of Hg remains predominantly atmospheric. This also compares well with existing records (Fig. 3). At Meishan, which was located in the margins of the Yangtze carbonate platform, the Hg and Hg/ TOC spike is coincident with more negative δ 202 Hg values (Fig. 1)  Hence, oxidation of terrestrial biomass is a compelling scenario to explain the palaeontological, sedimentological and geochemical data. There is clear observational evidence for the collapse of the terrestrial ecosystems and cessation of terrestrial C burial, stratigraphic evidence supporting the sequence and timing of the events (onset of the δ 13 C shift-collapse of the terrestrial ecosystem-Hg and C spike), sufficient quantity of Hg available, consistency with the isotopic evidence for changing Hg sources, and consistency with the δ 13 C records.
Massive terrestrial biomass oxidation during the PTME. Using our coupled C-Hg biogeochemical model, we show that the massive collapse of terrestrial ecosystems and oxidation of terrestrial biomass during the Permian-Triassic extinction had a huge impact on global Hg and C biogeochemistry. Hg stored in the terrestrial reservoirs was rapidly released as a consequence of the loss of terrestrial biomass and increased soil erosion 9,14 . This mechanism is the best explanation for the sharp increased loading of Hg into both terrestrial and marine water bodies and the negative shift in δ 202 Hg in coincidence with the marine mass extinction. Contemporaneously, increased soil carbon oxidation introduced large quantities of isotopically light C, accounting for the sharp negative δ 13 C anomaly registered in the sedimentary record (EP. I 30 ). In the model, the emission of Hg and C from magma and heating of sedimentary organic matter during the intrusive phase of the Siberian Traps LIP emplacement can account for the smaller,   two-to fourfold increase of Hg concentrations with respect to background levels, and the relatively longer negative δ 13 C trend that is recorded by both carbonates and organic matter, in marine and terrestrial settings.
A new scenario emerges for the PTME that links the collapse of ecosystems on land to the global geochemical changes recorded at the marine extinction interval. The disruption of terrestrial environments started during the initial phases of the Siberian Traps emplacement likely due to the release of volcanic gases as CO 2 , SO 2 and halogens, which could have triggered acid rain, ozone depletion, volcanic darkness, rapid cooling and subsequent global warming 8,49 . At the culmination of the terrestrial disturbance interval, when the ecosystems totally collapsed, large amounts of 13 C-depleted C and Hg deriving from a massive oxidation of terrestrial biomass were transported into aqueous habitats causing a steep decline in sedimentary δ 13 C (carbonates and organic matter), a sedimentary Hg concentrations spike and a shift in δ 202 Hg (Fig. 3). At this level, the marine mass extinction started. This, according to the existing chronostratigraphic framework, happened ∼60 Kyrs after the onset of the carbonisotope perturbation and of the terrestrial ecological disturbances 14 (Fig. 1).
The biogeochemical cycle of Hg is intimately linked to the cycle of organic matter and its constituting elements, such as C, N, S and P 50 . Hence, besides Hg and C, other organically-bound species would have been transferred from the terrestrial reservoirs into the marine system in large quantities at EP. I (Fig. 1). Addition of these species, particularly the nutrients P and N, are easily capable of driving ecosystem turnover, anoxia and eutrophication, and it is likely that this terrestrial input contributed to the marine extinction 9,11 . Our model does not include these additional cycles, but other models have shown that a relatively small increase in marine P delivery (2-3-fold) has the potential to drive marine anoxia or euxinia 51,52 . The scale of the terrestrial ecosystem collapse at the PTME could explain the severity of the biotic crisis at the Permian-Triassic boundary at all trophic levels, and should be a key consideration for future research. For other events, the Hg records are not so consistent nor as detailed as for the PTME. However, it is very likely that future research on other intervals could show the same Hg and C patterns as for the PTME.

Methods
Model derivation. This model is designed to track the transfer and isotopic signature of atmospheric and marine carbon and mercury over geological time, while being broadly applicable to changes on the timescale of ocean circulation. The biogeochemical system is taken largely from ref. 36 , with some additions from refs. 37,53,54 , with the underlying hydrological model from ref. 35 . The Hg cycle follows ref. 25 .
Model structure. The model has three ocean boxes: surface (s), high latitude (h) and deep (d). As in ref. 35 , the surface box is 100-m deep and occupies 85% of the ocean surface, whereas the high-latitude box is 250-m deep and represents 15% of the ocean surface. Each ocean box includes the same biogeochemical species, and a thermohaline circulation mixes the boxes in the order s, h, d. The upper boxes exchange with the atmosphere, which is a single box. As well as transfer fluxes between ocean and atmosphere boxes, biogeochemical fluxes of weathering, degassing and burial operate between the surface system and crust.
Model species. All model species are shown in Table 1.
Model fluxes. Model fluxes, with equations and present values are shown in Table 2.
Non-flux calculations. Atmospheric CO 2 volume ratio is calculated as: where CO 2a is atmospheric CO 2 in moles, and CO 2a 0 is this value at present day. Global average surface temperature (GAST) is: where k clim is climate sensitivity to doubling CO 2 , and t geol is time in millions of years before present and is expressed in negative terms. Low-latitude surface temperature (T s ) is assumed to scale by 2 3 times global temperature change, and both high-latitude (T h ) and deep (T d ) temperature are assumed to follow global temperature change.
For carbonate speciation, effective equilibrium constants are calculated following refs. 36,55 : K carb ¼ 5:75 10 À4 þ 6 10 À6 ðT j À 278Þ K CO 2 ¼ 0:035 þ 0:0019ðT j À 278Þ Dissolved carbon species are then calculated following Walker and Kasting 36 : We also explicitly calculate [H + ] concentration to observe model pH: ½CO 2À 3 Calcium carbonate saturation state is calculated as: where Ω j is the CaCO 3 saturation state in box j, and K sp is the solubility product. For terrestrial chemical weathering, temperature dependence of basalt and granite weathering is calculated as:   25 , and is assumed equivalent to the total Hg burial in their model to close the system over long timescales. This value is within their stated reasonable range but is~2× their chosen volcanic flux, the value is further increased by 50% due to higher rates of tectonic degassing through the Permian and Triassic. Runoff and burial fluxes are calculated to close the system under volcanic and wildfire input. Deposition and evasion fluxes directly follow ref. 25 . dt ¼ À f airsea s δ 13 C atm À f airsea h δ 13 C atm þ f ccdeg δ 13 C C þ f ocdeg δ 13 C G þ f oxidw δ 13 C G À f locb δ 13 C atm À ΔB À Á À f carbw δ 13 C atm À 2f silw δ 13 C atm þ f CO 2input δ 13 C input δ 13 C of low-latitude surface ocean DIC: d δ 13 DIC s Á DIC s À Á dt ¼ f airsea s δ 13 C atm þ tran DIC ds δ 13 DIC d À tran DIC sh δ 13 DIC s þ f carbw δ 13 C atm þ f carbw δ 13 C C þ 2f silw δ 13 C atm À f mccb δ 13 DIC s À f mocb ðδ 13 DIC s À ΔBÞ δ 13

Data availability
The geochemical data used in this paper come from already published literature, as cited in the text.