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.

Although the biodiversity goes up and down and the animal evolution has been frequently punctuated by mass extinction and the subsequent recovery 1 , the Earth System remained habitable for animals in the past 540 million years (Myr). The sustained habitability was warranted by the limited variations of the temperature and the redox condition of the surface Earth. It is widely accepted that the global temperature is mainly controlled by the pCO 2 level in the atmosphere 2 , while the surface Earth redox condition is determined by the atmosphere O 2 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 pO 2 level ever dropped below the threshold of 10% present atmospheric level (PAL) for animal breath 6 . It is proposed that the atmospheric O 2 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 pCO 2 and pO 2 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 pCO 2 and pO 2 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 phosphorus (P), 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 N 2 gas as ingredient 13 , the P mainly derives www.nature.com/scientificreports www.nature.com/scientificreports/ 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 CO 2 drawdown or O 2 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 seawater P would either be either converted to 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 (P in ) and the fraction of organophosphorus burial with respect to the total marine phosphorus burial (R p ) by using the seawater strontium isotope and the carbonate carbon isotope data.

Model Descriptions
Data acquisition. The Phanerozoic Strontium isotope (( ) ) Sr Sr sw 87 86 and carbonate carbon isotope (δ 13 C carb ) data were collected from chapters in the book: The Geologic Time Scale 2012 18,19 . Based on the variation of ( ) and δ 13 C carb , 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 water is the only major sources of both 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][21][22] . Overall, particulate P would have limited effects on the atmospheric pO 2 and pCO 2 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 CO 2 and consumes oxidants, offsetting, though not necessarily counterbalancing, the CO 2 consumption and O 2 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 (Fe 2+ ), resulting in the CO 2 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 pO 2 and pCO 2 levels. As such, primary productivity fueled by riverine dissolved P would directly impact the pCO 2 and pO 2 levels of the atmosphere.
To reconstruct the terrestrial dissolved P input, we approach with ( ) Sr Sr sw 87 86 , which documents the mixing between hydrothermal and riverine inputs 21,22 . Assuming the marine Sr budget (M Sr 0 ) 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: By solving the differential Eq. 1, riverine influx of Sr (F Sr r ) 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 (F P in ) and Sr inputs (F Sr r ) can be expressed by the following equation: where γ is the ratio between P and Sr concentrations of the weathered continental rocks, f dSr and f dp are the fraction of dissolved Sr and P fluxes in the total terrestrial Sr and P inputs.
Quantifying the marine phosphorus cycle. Seawater dissolved P can be buried with 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 org is the amount of organic carbon burial, R p P P org in = 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][24][25] .
The mass and carbon isotope balances can be expressed as: where Δ is the biological C isotope fractionation in photosynthesis, F C in and δ in are the flux and isotopic composition of DIC from continental weathering, F pre is the flux of marine carbonate precipitation, δ org is the isotopic composition of buried organic carbon, and δ DIC is the isotopic composition of the seawater DIC that is recorded in marine carbonate, i.e. δ δ = C C DIC c arb 13 13 . M 0 C is the size of dissolved inorganic carbon (DIC) pool in the ocean. The terrestrial DIC input is also related to the continental weathering, which is reflected by the mass ratio between the dissolution of 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 riverine DIC and dissolved Sr influxes, we define the following relationship: where A is the coefficient between mass of silicate rock weathering and riverine Sr input, and can be determined by the modern values (Table S2). Combining Eqs. 4-9, R p can be calculated by the following equation:   21 . Hydrothermal Sr influx (F Sr hy ) is estimated from the oceanic crust production rate (assuming the linear correlation), given the modern F Sr hy of 1 Tg/yr and the production rate of oceanic crust of 2.9 km 2 /yr 27,28 . δ CarW and δ SilW are assigned as −2.5‰ and −5‰ 29 . The ratio between F C in and F Sr r (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 carbonate rock weathering (F CarW ) N (130 Tg/yr, 30 2 ) 30 . The 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 × 10 18 g, i.e. double the DIC pool within the integrated interval, N (0,39 2 ) (10 18 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 km 2 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 simulations with repeated assigning the unconfined parameters (F CarW ΔM, and γ) 1000 times, and the mean value and standard deviation are calculated from 1000 iterations. The model calculate the modern F Sr r , P in and R p 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 (P in ) and fraction of organophosphorus burial (R p ) can be calculated from the Phanerozoic carbonate carbon and strontium isotope data (Fig. 1A,B, Table S1). The modeling results indicate that P in and R p display the long-term opposite trend (Fig. 1C,D). P in shows a general decreasing trend in Paleozoic, remains at a low value in Mesozoic, and increases in Cenozoic. In contrast, R P displays an increasing trend in Paleozoic followed by a high-value interval in Mesozoic, and declines in Cenozoic (Fig. 1C,D). Nevertheless, both P in and R P show more complex variations in the shorter time intervals. P in displays a general decreasing trend from ~2 Tg/yr in the early Cambrian to ~1 Tg/yr in late Ordovician (540-470 Ma), while R P remains at a low values of ~0.1 within this time interval. In the subsequent 40 million years (470-420 Ma), P in first increases to ~1.5 Tg/yr in late Ordovician, remains at the high value in early Silurian (445-425 Ma), and decreases to ~0.8 Tg/yr in the Silurian-Devonian boundary. In contrast, R p 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 P in from ~1 Tg/yr to ~0.5 Tg/yr, but R p shows more complex variations. R p 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), R p 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), P in slightly oscillates around ~0.5 Tg/yr, while R p shows a slight decrease from 0.3 to 0.2. P in increases to the modern value of 1.2 Tg/yr in the first 30 million years of Cenozoic (66-35 Ma), while R p displays an opposite trend. In the last 35 million years, both P in and R P remains at the modern values (Fig. 1C,D). www.nature.com/scientificreports www.nature.com/scientificreports/ coupling of terrestrial input of dissolved p and marine p cycle. The negative correlation between P in and R p implies the coupling of terrestrial P input and marine P cycle (Figs. 1C,D, 2A). An increase of P in would associate with the lower efficiency of biological P utilization in the ocean ( Fig. 2A). In contrast, when the riverine influx of dissolved P is low, more P is buried as organophosphorus, implying less fraction of dissolved P being sequestered as authigenic phosphate or iron-bound P 15 .
The riverine input of dissolved P determines the upper bound of marine primary productivity, and confines the maximum burial of 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 iron (Fe) redox cycle plays the key role in the P biogeochemical cycle between seawater and sediment 15,17 . In oxic conditions, ferrous Fe (Fe 2+ ) is thermodynamically unstable and would be spontaneously oxidized to ferric Fe (Fe 3+ ). 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 Ca 2+ and HCO 3 − in sediment porewater 17,36 . On the other hand, Fe 2+ 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 Fe redox cycle is mainly driven by the microbial iron reduction (MIR) that uses reactive ferric Fe (e.g. Fe 2 O 3 , 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 O 2 fugacity at the sediment-water interface (SWI), enhancing the benthic Fe 2+ flux from sediment porewater to seawater [37][38][39][40] . Therefore, high P in favors repeated Fe reduction-oxidation cycles between seawater and sediment, enhancing the shuttling of seawater P into sediment.
The negative correlation between P in and R P also implies the negative feedbacks imposed by marine P cycle as a response to the change of riverine input of dissolved P (Fig. 1), which limits the variation of organic carbon burial. As shown in Fig. 1, though δ 13 C carb displays dramatic oscillation in Phanerozoic (Fig. 1A), the amount of organic carbon burial (F org ) remains more or less constant throughout Phanerozoic (Fig. 1E).
Although generally negatively correlated, P in and R p do not show exact opposite trends (Fig. 1C,D). For example, in late Ordovician and Silurian, P in 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, R p displays a steady increase from 0.08 to 0.25 (Fig. 1C,D). The reason for the sudden increase of P in 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 R p 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][43][44] , and (2) the ocean oxygenation in Silurian as evidenced by the global formation of marine red beds depositions 45,46 . www.nature.com/scientificreports www.nature.com/scientificreports/ In addition, a decrease of P in in late Cretaceous (100-65 Ma) does not associate with an increase of R P . Instead, R P remains nearly invariant throughout Cretaceous. The invariant R P might be attributed to, in the context of super-greenhouse climatic condition in Cretaceous, the widespread oceanic anoxia 47,48 , under which condition the 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 indicates a decrease of P in in late Paleozoic (Fig. 1C), when vascular 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 CO 2 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 CO 2 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 riverine influx of P not only depends on the rate of chemical weathering, but also is a function of physical 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 P in 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 P in 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 P in 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 .
Diversification of vascular land plant in late Devonian to early Carboniferous (400-340 Ma) is also witnessed the significant decoupling of P in and R P . R p displays a sharp increase from ~0.1 to 0.3 followed by a decrease to 0.08, while P in shows a slight decrease from 1.2 to 0.8 Tg/yr (Fig. 1C,D). The reason for the rapid change of R P is unknown. We speculate that the increase of R p 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 P in reflects the long-term change of continental weathering, while R p is driven by the short-term marine P cycles. Organic carbon burial as well as δ 13 C carb was controlled by both the long-term and short-term processes. Therefore, δ 13 C carb does not show any correlation with either P in or R p (Fig. 1A,C,D). In addition, there is no correlation between F org and P in (yellow circles in Fig. 2B) or between F org and R p (yellow circles in Fig. 2C) for the long duration intervals (>3 Myr). In contrast, the short-term (<3 Myr) δ 13 C carb excursions might be mainly controlled by the change of R p . This argument is well supported by the crude correlation between F org and R p in the short duration intervals (black dots in Fig. 2C).
Some short duration δ 13 C carb excursions, such as, the negative δ 13 C carb excursion 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 R p value (Table S1). The negative R p 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 indicates that the continent and ocean systems have been tightly coupled throughout the Phanerozoic. Though the continent system could be perturbed by, e.g., assembly and breakup of supercontinents 56 , eruption of large igneous provinces 61 , or evolution of terrestrial ecosystem 10 , resulting in a dramatic change of 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 pCO 2 and pO 2 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.

Data availability
All data is available in the main text or the supplementary materials.