The coupling of Phanerozoic continental weathering and marine phosphorus cycle.

Organic matter production and decomposition primarily modulate the atmospheric O2 and CO2 levels. The long term marine primary productivity is controlled by the terrestrial input of phosphorus (P), while the marine P cycle would also affect organic matter production. In the past 540 million years, the evolution of terrestrial system, e.g. colonization of continents by vascular land plants in late Paleozoic, would certainly affect terrestrial P input into the ocean, which in turn might have impacted the marine primary productivity and organic carbon burial. However, it remains unclear how the marine P cycle would respond to the change of terrestrial system. Here we reconstruct the secular variations of terrestrial P input and biological utilization of seawater P in Phanerozoic. Our study indicates that riverine dissolved P input and marine P biological utilization (i.e. the fraction of P being buried as organophosphorus) are inversely correlated, suggesting the coupling of continental P input and marine P cycle. We propose an increase of P input would elevate surface ocean productivity, which in turn enhances marine iron redox cycle. Active Fe redox cycle favors the scavenging of seawater P through FeOOH absorption and authigenic phosphate formation in sediments, and accordingly reduces the bioavailability of seawater P. The negative feedback of marine P cycle to terrestrial P input would keep a relatively constant organic carbon burial, limiting the variations of surface Earth temperature and atmospheric O2 level.

Introduction

Although the biodiverclass="Chemical">sity goes uclass="Chemical">p and down and the animal evolution has been frequently class="Chemical">punctuated by mass extinction and the subsequent recovery[1], the Earth System remained habitable for animals in the class="Chemical">past 540 million years (Myr). The sustained habitability was warranted by the limited variations of the temclass="Chemical">perature and the class="Chemical">pan class="Disease">redox condition of the surface Earth. It is widely accepted that the global temperature is mainly controlled by the pCO2 level in the atmosphere[2], while the surface Earth redox condition is determined by the atmosphere O2 content[3]. In spite of the long-term oscillation between the icehouse and greenhouse/hothouse climatic conditions[4], devastating extreme climate events, such as the global glaciations, have never occurred after the Cryogenian (~720-635 Ma) snowball Earth events[5]. Nor has the atmospheric pO2 level ever dropped below the threshold of 10% present atmospheric level (PAL) for animal breath[6]. It is proposed that the atmospheric O2 level might have increased from 10–40% PAL in the early Cambrian to the modern level in late Paleozoic, and showed limited variation in the past 250 million years[6,7], although oceanic anoxia could be common and widespread in certain geological ages[8].

The atmospheric class="Chemical">pCO2 and class="Chemical">pan class="Chemical">pO2 levels are modulated by the global carbon cycle, in which organic matter production and destruction in both terrestrial and marine systems are the core processes[9]. The terrestrial carbon cycle might be negligible before the diversification of vascular land plant in late Paleozoic[10,11], but the marine carbon cycle has been the key component of global carbon cycle throughout the Earth’s history. As such, the secular variation of marine primary productivity would provide a key constraint on the evolution of atmospheric pCO2 and pO2 levels.

The marine primary productivity is dominated in the upper most 200 m of the ocean (the euphotic zone), and is controlled by the availability of macronutrients, including class="Chemical">phosphorus (P), class="Chemical">pan class="Chemical">nitrogen (N) and silicon (Si), and micronutrients, such as iron (Fe), zinc (Zn) etc.[12]. Unlike the bio-available N that can be synthesized by diazotrophs (nitrogen-fixing cyanobacteria) by using atmospheric N2 gas as ingredient[13], the P mainly derives from terrestrial input[14,15]. It has long been held that the riverine P input primarily controls the marine primary productivity in the geological time scale[15,16], placing the upper bound of surface ocean primary productivity and confining the largest magnitudes of CO2 drawdown or O2 rise in the atmosphere.

On the other hand, marine P cycle might impact the short-term primary productivity at local or regional scales[15,16]. Substantial amount of seaclass="Chemical">water P would either be either converted to class="Chemical">pan class="Chemical">authigenic phosphate or absorbed by iron oxyhydroxides (FeOOH)[15,17], reducing the bio-availability of seawater P. In the modern abyssal sediments, authigenic phosphate precipitation and absorption by FeOOH would account for more than 80% of P burial[15]. Thus, the long-term marine primary productivity is controlled by both riverine P influx and marine P cycle. In this study, we reconstruct the Phanerozoic secular variations of terrestrial input of dissolved P (Pin) and the fraction of organophosphorus burial with respect to the total marine phosphorus burial (Rp) by using the seawater strontium isotope and the carbonate carbon isotope data.

Model Descriptions

Data acquisition

The Phanerozoic class="Chemical">Strontium isotoclass="Chemical">pe and class="Chemical">pan class="Chemical">carbonate carbon isotope (δ13Ccarb) data were collected from chapters in the book: The Geologic Time Scale 2012[18,19]. Based on the variation of and δ13Ccarb curves, Phanerozoic (541 Ma-present) was divided into 65 intervals, each of which is defined by distinct trends of and/or δ13Ccarb as compared with adjacent intervals (Table S1). The time intervals can be classified as the long (>3 Myr) and short (= or <3 Myr) duration intervals, reflecting the long-term and short-term variations in and δ13Ccarb, respectively. Totally 41 long duration intervals and 24 short duration intervals were identified (Table S1), and the total range of the long duration intervals (493 Myr) accounts for 91% of Phanerozoic (541 Myr).

Quantifying the riverine influx of dissolved phosphorus

River class="Chemical">water is the only major sources of both class="Chemical">pan class="Chemical">Sr and P in the ocean, and the riverine influxes are directly controlled by the continental weathering. Riverine P input includes dissolved and particulate P, while Sr is mainly transported in its dissolved form[15,18,20-22]. Overall, particulate P would have limited effects on the atmospheric pO2 and pCO2 levels for the following reasons. Firstly, detrital particulate inorganic P (e.g. unweathered apatite from igneous rocks or fossil authigenic phosphate from sedimentary rocks) is geochemically inactive, and thus cannot involve in the biogeochemical reactions in the ocean. Secondly, although particulate organic P could be mobilized by the degradation of associated organic component and might have different C: P ratios from the Redfield ratio of 106: 1[23], organic matter remineralization releases CO2 and consumes oxidants, offsetting, though not necessarily counterbalancing, the CO2 consumption and O2 emission in photosynthesis. Thirdly, the iron-bound P from continental input could be remobilized by reduction of ferric Fe, which is associated with organic matter remineralization and releases of ferrous Fe (Fe2+), resulting in the CO2 emission and the reduction of oxidative state of the ocean[14,15,20]. Thus, remobilization of particulate P tends to cancel the effect of organic matter production and thus have limited impact on the atmospheric pO2 and pCO2 levels. As such, primary productivity fueled by riverine dissolved P would directly impact the pCO2 and pO2 levels of the atmosphere.

To reconstruct the terrestrial dissolved P input, we approach with , which documents the mixing between hydrothermal and class="Chemical">riverine inclass="Chemical">puts[21,22]. Assuming the marine class="Chemical">pan class="Chemical">Sr budget () remains a constant, i.e. the riverine and hydrothermal Sr inputs are counterbalanced by carbonate precipitation, the isotope mass balance of marine Sr cycle can be expressed as:where is the Sr flux of source or sink i, and the subscripts r and hy represent the riverine and hydrothermal inputs, respectively. Because carbonate precipitation is the only major sink of seawater Sr, and Sr in carbonate minerals records the seawater Sr isotopic composition, i.e. [21]

By solving the difclass="Chemical">ferential Eq. 1, class="Chemical">pan class="Chemical">riverine influx of Sr () can be calculated by the following equation:where subscript sw1 and sw2 are the seawater compositions at the beginning and the end of the time interval. Because the riverine dissolved P is mainly derived from weathering of continent rocks[15], the relationship between the riverine dissolved P () and Sr inputs () can be expressed by the following equation:where γ is the ratio between P and Sr concentrations of the weathered continental rocks, f and f are the fraction of dissolved Sr and P fluxes in the total terrestrial Sr and P inputs.

Quantifying the marine phosphorus cycle

Seaclass="Chemical">water dissolved P can be buried with class="Chemical">pan class="Chemical">organic carbon as organophosphorus, as authigenic phosphate, and by adsorption onto FeOOH, i.e. Fe-bound phosphorus[14,15,17]. The amount of organophosphorus burial is proportional to the quantity of organic matter burial, which can be represented by the following equation:where F is the amount of organic carbon burial, represents the proportion of organophosphorus burial with respect to the total dissolved P input from continental weathering, and RED is the Redfield ratio of marine primary productivity with a theoretical C: P (molar ratio) of 106:1[23-25].

The mass and pan class="Chemical">carbon isotoclass="Chemical">pe balances can be exclass="Chemical">pressed as:

where Δ is the biological C isotope fractionation in photosyntheclass="Chemical">sis, and are the flux and isotoclass="Chemical">pic comclass="Chemical">poclass="Chemical">pan class="Chemical">sition of DIC from continental weathering, is the flux of marine carbonate precipitation, is the isotopic composition of buried organic carbon, and is the isotopic composition of the seawater DIC that is recorded in marine carbonate, i.e. . M0C is the size of dissolved inorganic carbon (DIC) pool in the ocean.

The terrestrial class="Chemical">DIC inclass="Chemical">put is also related to the continental weathering, which is reflected by the mass ratio between the dissolution of class="Chemical">pan class="Chemical">silicate rock and carbonate rock[26]. The mass and isotopic balance of continent weathering can be expressed as:where the subscript CarW and SilW represent carbonate rock and silicate rock weathering, respectively.

To link the class="Chemical">riverine class="Chemical">pan class="Chemical">DIC and dissolved Sr influxes, we define the following relationship:

where A is the coefficient between mass of class="Chemical">silicate rock weathering and class="Chemical">pan class="Chemical">riverine Sr input, and can be determined by the modern values (Table S2). Combining Eqs. 4–9, R can be calculated by the following equation:

Parameter setting

The Phanerozoic seaclass="Chemical">water class="Chemical">pan class="Chemical">Sr isotope and carbonate carbon isotope data are tabulated in Table S1, and the assigned values of all parameters are listed in Table S2. , , and are 0.7119, 0.7035, and 1.12 × 1019 g, respectively[21]. Hydrothermal Sr influx () is estimated from the oceanic crust production rate (assuming the linear correlation), given the modern of 1 Tg/yr and the production rate of oceanic crust of 2.9 km2/yr[27,28]. and are assigned as −2.5‰ and −5‰[29]. The ratio between and (A) is calculated from the modern values of 89[30]. The biological fractionation of organic carbon production (Δ) is 25‰[31].

We further assume the normal distribution of the class="Chemical">carbonate rock weathering () N (130 Tg/yr, 302)[30]. The class="Chemical">pan class="Chemical">size of seawater DIC could also fluctuation. We assume that normal distribution of the rate of seawater DIC pool variation (△M) with the sigma of 39 × 1018 g, i.e. double the DIC pool within the integrated interval, N (0,392) (1018 g). The parameter γ, the P/Sr ratio of weathered continental rocks is dependent on the lithology of exposed rocks. Especially, the exposure of large igneous province (LIP) in continents, e.g. the Siberian trap with the area of 2.1 million km2 [32], would el evate the P/Sr ratio in certain geological interval. We assume the probability of LIP eruption with the exposure area of Siberian trap is 0.01, and the normal dis tribution of γ with the mean value of the upper continental crust composition (P: 655 ppm, Sr: 320 ppm)[33].

We run Monte Caro class="Chemical">simulations with reclass="Chemical">peated asclass="Chemical">pan class="Chemical">signing the unconfined parameters (F ΔM, and γ) 1000 times, and the mean value and standard deviation are calculated from 1000 iterations. The model calculate the modern , Pin and Rp of 2.9 Tg/yr, 1.2 Tg/yr and 0.15, respectively, consistent with the observed values[15,21].

Secular variation of terrestrial input of dissolved P and marine P cycle

The secular variations of terrestrial input of dissolved P (Pin) and fraction of class="Chemical">organophosphorus burial (Rclass="Chemical">p) can be calculated from the class="Chemical">pan class="Chemical">Phanerozoic carbonate carbon and strontium isotope data (Fig. 1A,B, Table S1). The modeling results indicate that Pin and Rp display the long-term opposite trend (Fig. 1C,D). Pin shows a general decreasing trend in Paleozoic, remains at a low value in Mesozoic, and increases in Cenozoic. In contrast, RP displays an increasing trend in Paleozoic followed by a high-value interval in Mesozoic, and declines in Cenozoic (Fig. 1C,D).

Figure 1

(A) Secular variation of carbonate carbon isotopes. Data are from Saltzman and Thomas[19] and Prokoph, et al.[64]. (B) Secular variation of seawater strontium isotopes. Data are from McArthur, et al.[18]. (C) Secular variation of riverine influx of dissolved phosphate (Pin). (D) Secular variation of fraction of organophosphorus burial with respect to total riverine influx of dissolved P (Rp). (E) Secular variation of mass of organic carbon burial (Forg). (F) Evolutionary trend of silica-secreting organisms (diatom and radiolarian), showing the stepwise diversification in Mesozoic[62], and the width of the box represents the diversity; also marked the time lines of vascular land plants evolution[11] and assembly and breakup of Pangea[63]. The solid lines in C,D,E are the mean values of Monte Caro simulation with 1000 iterations, while the shadowed regions bracket the 2 times of standard deviations (i.e. 95% confidence interval) of 1000 iterations.

(A) Secular variation of class="Chemical">carbonate carbon isotoclass="Chemical">pes. Data are from Saltzman and Thomas[19] and Prokoclass="Chemical">ph, et al.[64]. (B) Secular variation of seaclass="Chemical">pan class="Chemical">water strontium isotopes. Data are from McArthur, et al.[18]. (C) Secular variation of riverine influx of dissolved phosphate (Pin). (D) Secular variation of fraction of organophosphorus burial with respect to total riverine influx of dissolved P (Rp). (E) Secular variation of mass of organic carbon burial (Forg). (F) Evolutionary trend of silica-secreting organisms (diatom and radiolarian), showing the stepwise diversification in Mesozoic[62], and the width of the box represents the diversity; also marked the time lines of vascular land plants evolution[11] and assembly and breakup of Pangea[63]. The solid lines in C,D,E are the mean values of Monte Caro simulation with 1000 iterations, while the shadowed regions bracket the 2 times of standard deviations (i.e. 95% confidence interval) of 1000 iterations.

Nevertheless, both Pin and RP show more complex variations in the shorter time intervals. Pin displays a general decreaclass="Chemical">sing trend from ~2 Tg/yr in the early Cambrian to ~1 Tg/yr in late Ordovician (540-470 Ma), while RP remains at a low values of ~0.1 within this time interval. In the subsequent 40 million years (470-420 Ma), Pin first increases to ~1.5 Tg/yr in late Ordovician, remains at the high value in early class="Chemical">pan class="Chemical">Silurian (445-425 Ma), and decreases to ~0.8 Tg/yr in the Silurian-Devonian boundary. In contrast, Rp illustrates a steady increase from 0.1 to 0.2 within 40 million years (470-430 Ma) and a sharp decline to 0.1 in the following 30 million years (430-400 Ma). In the following 150 million years (400-250 Ma), there is a continuous decline of Pin from ~1 Tg/yr to ~0.5 Tg/yr, but Rp shows more complex variations. Rp increases sharply from 0.1 to 0.3 from 400 Ma to 370 Ma and then decreases to 0.1 in the following 30 million years (370-340 Ma). In the last 40 million years of Carboniferous (340-300 Ma), Rp returns to the highest value of 0.3 and remained in the plateau throughout Permian until the Permian-Triassic transition that is characterized by a rapid decline to 0.2. In the whole Mesozoic (250-65 Ma), Pin slightly oscillates around ~0.5 Tg/yr, while Rp shows a slight decrease from 0.3 to 0.2. Pin increases to the modern value of 1.2 Tg/yr in the first 30 million years of Cenozoic (66-35 Ma), while Rp displays an opposite trend. In the last 35 million years, both Pin and RP remains at the modern values (Fig. 1C,D).

Coupling of terrestrial input of dissolved P and marine P cycle

The negative correlation between Pin and Rp implies the coupling of terrestrial P input and marine P cycle (Figs. 1C,D, 2A). An increase of Pin would associate with the lower efficiency of biological P utilization in the ocean (Fig. 2A). In contrast, when the class="Chemical">riverine influx of dissolved P is low, more P is buried as class="Chemical">pan class="Chemical">organophosphorus, implying less fraction of dissolved P being sequestered as authigenic phosphate or iron-bound P[15].

Figure 2

(A) Crossplot showing the negative correlation between Pin and Rp. (B) Crossplot showing the absence of correlation between Forg and Pin. (C) Crossplot showing the absence of correlation between Forg and Rp for the long duration intervals (yellow circles), and a crude correlation for the short duration intervals (dark dots). The yellow circles represent the intervals greater than 3 million years (long duration) while the dark dots indicate intervals less than 3 million years (short duration).

(A) Crossplot showing the negative correlation between Pin and Rp. (B) Crossplot showing the absence of correlation between class="Chemical">Forg and Pin. (C) Crossclass="Chemical">plot showing the absence of correlation between class="Chemical">pan class="Chemical">Forg and Rp for the long duration intervals (yellow circles), and a crude correlation for the short duration intervals (dark dots). The yellow circles represent the intervals greater than 3 million years (long duration) while the dark dots indicate intervals less than 3 million years (short duration).

The class="Chemical">riverine inclass="Chemical">put of dissolved P determines the uclass="Chemical">pclass="Chemical">per bound of marine class="Chemical">primary class="Chemical">productivity, and confines the maximum burial of class="Chemical">pan class="Chemical">organic carbon. However, not all seawater P would be buried with organic matter, instead, substantial amount of seawater P would be eventually sequestered as authigenic phosphate or iron-bound P[15,17], reducing the bioavailability of seawater P. In the modern abyssal sediment where detrital P deposition is negligible, authigenic phosphate and iron-bound P account for ~70% and ~15% of P burial, respectively[15].

The class="Chemical">iron (class="Chemical">pan class="Chemical">Fe) redox cycle plays the key role in the P biogeochemical cycle between seawater and sediment[15,17]. In oxic conditions, ferrous Fe (Fe2+) is thermodynamically unstable and would be spontaneously oxidized to ferric Fe (Fe3+). In the oxic seawater with near neutral pH of 8–8.2, ferric Fe would precipitate as iron oxyhydroxides (FeOOH), during which seawater P would be absorbed by FeOOH. The degree of P absorption, measured by the P to Fe molar ratio [(P/Fe)100] of iron oxides/oxyhydroxides deposits, is controlled by the concentrations of dissolved P and Si in seawater, the latter of which would compete with P for the absorption sites in FeOOH[34,35]. Precipitation of FeOOH particles shuttles seawater P into sediments. Further reduction of FeOOH by iron reducing microbes (IRM) in sediment liberates absorbed P into porewater, resulting in the precipitation of authigenic phosphate by reacting with Ca2+ and HCO3− in sediment porewater[17,36]. On the other hand, Fe2+ could diffuse back into seawater, i.e. the benthic iron flux[37,38], where ferrous Fe oxidation and FeOOH precipitation would continue shuttling seawater P into sediments[34]. If FeOOH is not reduced in sediment, e.g. due to the depletion of organic matter, preservation of FeOOH would result in the Fe-bound P burial. Therefore, the amount of iron-bound P burial is determined by the availability of reactive Fe that involves in the redox cycle between seawater and sediment, while authigenic phosphate burial is controlled by both the number of Fe redox cycling and the amount of cycled reactive Fe.

The class="Chemical">Fe class="Chemical">pan class="Disease">redox cycle is mainly driven by the microbial iron reduction (MIR) that uses reactive ferric Fe (e.g. Fe2O3, FeOOH) as electron acceptor to oxidize organic matter to bicarbonate[36]. Thus, active Fe redox cycle requires sufficient supplies of organic matter. The Fe redox cycle in the ocean and continental weathering are related in that, high riverine flux of dissolved P would stimulate high surface ocean productivity[15], which in turn enhances organic matter supply for MIR. On the other hand, high productivity also lowers O2 fugacity at the sediment-water interface (SWI), enhancing the benthic Fe2+ flux from sediment porewater to seawater[37-40]. Therefore, high Pin favors repeated Fe reduction-oxidation cycles between seawater and sediment, enhancing the shuttling of seawater P into sediment.

The negative correlation between Pin and RP also implies the negative class="Chemical">feedbacks imclass="Chemical">posed by marine P cycle as a resclass="Chemical">ponse to the change of class="Chemical">pan class="Chemical">riverine input of dissolved P (Fig. 1), which limits the variation of organic carbon burial. As shown in Fig. 1, though δ13Ccarb displays dramatic oscillation in Phanerozoic (Fig. 1A), the amount of organic carbon burial (Forg) remains more or less constant throughout Phanerozoic (Fig. 1E).

Although generally negatively correlated, Pin and Rp do not show exact oppoclass="Chemical">site trends (Fig. 1C,D). For examclass="Chemical">ple, in late Ordovician and class="Chemical">pan class="Chemical">Silurian, Pin first decreases from 1.5 to 1 Tg/yr, then followed by a sharp increase to  2 Tg/yr and decrease within 40 million years. In contrast, Rp displays a steady increase from 0.08 to 0.25(Fig. 1C,D). The reason for the sudden increase of Pin might be related to the closure of Iapetus ocean and the formation of Caledonian Orogeny[41], the latter of which might have substantially enhanced weathering input of P, while the steady increase of Rp might reflect, (1) the diversifications of phytoplanktons and zooplanktons in the Great Ordovician Biodiversification Event, which might have enhanced the efficiency of particulate organic carbon burial[42-44], and (2) the ocean oxygenation in Silurian as evidenced by the global formation of marine red beds depositions[45,46].

In addition, a decrease of Pin in late Cretaceous (100-65 Ma) does not associate with an increase of RP. Instead, RP remains nearly invariant throughout Cretaceous. The invariant RP might be attributed to, in the context of super-greenhouse climatic condition in Cretaceous, the widespread class="Disease">oceanic anoxia[47,48], under which condition the class="Chemical">pan class="Chemical">Fe redox cycle might be insensitive to the surface ocean productivity. It should be also noted that the flourish of silica-secreting organisms, such as diatom and radiolarian (Fig. 1F), might have lowered the seawater Si concentration[49], enhanced P absorption onto FeOOH particles[35], and accordingly promoted the efficiency of P delivery by FeOOH shuttle. However, it remains unclear how change of seawater Si concentration could affect the marine P cycle in Cretaceous.

Effect of vascular land plant diversification on the riverine P input and marine P cycle

The modeling result inclass="Chemical">dicates a decrease of Pin in late Paleozoic (Fig. 1C), when vascular class="Chemical">pan class="Species">land plants diversified and the forest began to form[10,11]. It is widely accepted that the biological weathering is more efficient than inorganic chemical weathering that uses CO2 as reactant[10,50]. With the presence of land plant, weathering intensity could be substantially enhanced in several ways. First of all, the decomposition of organic matter in soils would generate more corrosive organic acid that accelerates the reaction rate, while the respiration of plant roots also elevates CO2 concentration in soils, further enhancing mineral dissolution[51]. Moreover, plant roots would split rocks apart, favoring water percolation and increasing the reactive surface within rock interior. The enhanced continental weathering intensity in late Paleozoic is supported by the widespread sedimentary bauxite depositions since late Devonian[52]. Therefore, it is predicted that riverine influx of P would be higher, when continents were colonized by vascular land plants.

In fact, the class="Chemical">riverine influx of P not only declass="Chemical">pends on the rate of chemical weathering, but also is a function of class="Chemical">phyclass="Chemical">pan class="Chemical">sical erosion[53,54]. Rapid erosion would reduce the reaction time of chemical weathering. The erosion rate would be smaller in well-vegetated continents. and is sensitive to regional/local precipitation and temperature. We speculate that low Pin in late Paleozoic might be related to the amalgamation of supercontinents Pangea (Fig. 1C,F), during which the weakened hydrological cycles in the ocean-continent system would result in an expansion of drought area in the continent interior[55,56]. Low Pin in Carboniferous and Permian might also relate to the global cooling during the late Paleozoic ice age[57], since weathering rate is temperature dependent[58]. For the same reason, an increase of Pin in late Mesozoic might be attributed to an elevated continental hydrological cycle during the breakup of Pangea (Fig. 1C,F)[56] as well as the greenhouse/hothouse climatic conditions[4].

Diverclass="Chemical">sification of vascular land class="Chemical">plant in late Devonian to early class="Chemical">pan class="Chemical">Carboniferous (400-340 Ma) is also witnessed the significant decoupling of Pin and RP. Rp displays a sharp increase from ~0.1 to 0.3 followed by a decrease to 0.08, while Pin shows a slight decrease from 1.2 to 0.8 Tg/yr (Fig. 1C,D). The reason for the rapid change of RP is unknown. We speculate that the increase of Rp in mid-Devonian might be attributed to, e.g. oxygenation in atmosphere due to the formation of forest[6,11].

The long-term and short-term coupling of continent-ocean system

We propose that Pin reflects the long-term change of continental weathering, while Rp is driven by the short-term marine P cycles. class="Chemical">Organic carbon burial as well as δ13Ccarb was controlled by both the long-term and short-term class="Chemical">processes. Therefore, δ13Ccarb does not show any correlation with either Pin or Rclass="Chemical">p (Fig. 1A,C,D). In addition, there is no correlation between class="Chemical">pan class="Chemical">Forg and Pin (yellow circles in Fig. 2B) or between Forg and Rp (yellow circles in Fig. 2C) for the long duration intervals (>3 Myr). In contrast, the short-term (<3 Myr) δ13Ccarb excursions might be mainly controlled by the change of Rp. This argument is well supported by the crude correlation between Forg and Rp in the short duration intervals (black dots in Fig. 2C).

Some short duration δ13Ccarb excurclass="Chemical">sions, such as, the negative δ13Ccarb excurclass="Chemical">pan class="Chemical">sion between 525-523 Ma (the DICE event)[59], cannot be explained by the canonical carbon cycle model[29]. The DICE event is characterized by a negative Rp value (Table S1). The negative Rp may indicate additional P supply from a non-terrestrial source, such as remineralization of a dissolved organic carbon (DOC) pool in the deep ocean[60].

Closing remarks

Our study inclass="Chemical">dicates that the continent and ocean systems have been tightly couclass="Chemical">pled throughout the Phanerozoic. Though the continent system could be class="Chemical">perturbed by, e.g., assembly and breakuclass="Chemical">p of suclass="Chemical">percontinents[56], eruclass="Chemical">ption of large igneous class="Chemical">provinces[61], or evolution of terrestrial ecosystem[10], resulting in a dramatic change of class="Chemical">pan class="Chemical">riverine input of dissolved P. Such perturbations could be mitigated by marine P cycle. More specifically, the coupled continent and ocean systems would limit the variation of organic carbon burial, which plays the key role in modulating atmospheric pCO2 and pO2 levels[29], and accordingly the global temperature and the Earth's surface redox conditions. Thus, the coupling of continental P weathering and marine P cycle warrants the long-term stability and habitability of the Earth system in the past 540 million years, paving the way for the evolution of intelligent human beings.

Supplementary information.