Soil organic phosphorus transformations during 2000 years of paddy-rice and non-paddy management in the Yangtze River Delta, China

The contents and properties of soil organic phosphorus (Po) largely drive ecosystem productivity with increasing development of natural soil. We hypothesized that soil Po would initially increase with paddy management and then would persist under steady-state conditions. We analyzed soils from a 2000-year chronosequence of a rice-wheat rotation and an adjacent non-paddy 700-year chronosequence in Bay of Hangzhou (China) for their Po composition using solution 31P-NMR after NaOH-EDTA extraction. Land reclamation promoted Po accumulation in both paddy and non-paddy topsoils (depths ≤ 18 cm) until steady-state equilibria were reached within 200 years of land use. Greater Po concentrations were found, however, in the non-paddy subsoils than in those under paddy management. Apparently, the formation of a dense paddy plough pan hindered long-term Po accumulation in the paddy subsoil. The surface soils showed higher proportions of orthophosphate diesters under paddy than under non-paddy management, likely reflecting suppressed decomposition of crop residues despite elevated microbial P compounds stocks under anaerobic paddy-rice management. Intriguingly, the composition of Po was remarkably stable after 194-years of paddy management and 144-years of non-paddy management, suggesting novel steady-state equilibria of P dynamics had been reached in these man-made ecosystems after less than two centuries.

only 700-to 1000-years, the soils are decalcified and have significant iron (Fe) oxide transformations, which can alter the sorption of specific P o forms 14,16 .
Soil P o is comprised of a variety of compounds that differ in their stability and biological availability in environment 17 . The variation of P composition over pedogenic time scales has ecological significance 5 , as it may alter species composition of the seedling community, change the competitive ability of adult plant species 18 , and modulate soil microbial P cycling 19,20 . To reveal changes in P o composition, solution 31 P-nuclear magnetic resonance ( 31 P-NMR) spectroscopy after NaOH-Na 2 EDTA extraction has been the method of choice [21][22][23] . This allows the identification of specific P o species, such as deoxyribonucleic acid (DNA), αand β-glycerophosphate, mononucleotides, phosphonates, and even inorganic long-chain polyphosphate and pyrophosphate structures, including compounds originating from living microbial cells and fungal tissue 5,[24][25][26] . However, little information is available on the degree and rates at which different P o species reflect impacts of prolonged paddy or non-paddy management on the overall P o status of the soils.
The objective of this study was to evaluate how prolonged paddy management influenced the degree and rates at which various P o compounds accumulated in soil. For this purpose, we extracted P o with NaOH-Na 2 EDTA and characterized compounds by solution 31 P-NMR spectroscopy. Samples originated from a unique 2000-year-old chronosequence of paddy soil management; adjacent sites under non-paddy management for 700 years served as a reference. All soils had developed from tidal wetland sediments in the Yangtze River delta, China. We hypothesized that paddy management accelerated the accumulation of P o compounds in soil, which may then persist for several hundred years under steady-state conditions.

Results
During paddy soil management, the puddling of surface soil accompanies the formation of a dense plough (Ardp) horizon, which potentially restricts root growth into the deeper subsoil 27 . With increasing time of paddy soil development, as these Ardp horizons mature, topsoil surfaces may become increasingly decoupled from subsoil biogeochemistry 13,14 . In order to understand the role of paddy management in the biological cycling of P, particularly when compared to native ecosystems, we thus have to understand both the transformation of P with increasing soil depth as well as with time of management.
Changes in phosphorus speciation with soil depth. Total P concentrations decreased with increasing depth from ~800 mg P kg −1 in the 100-year-old paddy topsoils and from ~1100 mg P kg −1 in the respective non-paddy topsoils down to 40-50 cm soil depth (600-690 mg P kg −1 ) and then remained relatively constant in the subsoil ( Table 1). The concentrations of NaOH-Na 2 EDTA extractable total P (E P ) and inorganic P (E P -P i ) in these soils showed a similar trend ( Table 1). The contributions of extracted P o (E p -P o ) to total P in the topsoil of the paddy soil were almost double those under non-paddy management, whereas those in paddy subsoil were reduced compared with subsoil under non-paddy management ( Table 1). The contribution of E p -P o declined with depth while that of E p -P i increased ( Table 1). As a result, all P o species in the paddy subsoil (from 30 to 100 cm) showed very low 31 P-NMR signal intensities, and the spectra were dominated by orthophosphate P ( Table 2). The non-paddy site showed a greater proportion of total P as orthophosphate monoesters below 30 cm soil depth than the paddy site, but the lowest subsoil was again dominated by P i (Table 2).
Above 30 cm depth, orthophosphate diesters and myo-IHP were much higher proportions of E P in the paddy soils than in the non-paddy soils. There was a substantial increase in phosphonates at the 30-50 cm depth in the paddy soils. Overall, the effects of paddy management on P o dynamics were thus mainly detectable for the surface soils. Hence, we hereafter concentrated on topsoil samples (~0-15 cm) in our study of P transformations along the chronosequence.  Table 1. Bulk density and concentrations of total P (P-total), NaOH-Na 2 EDTA extracted total P (E p ), organic carbon (OC) and total nitrogen (N-total) at different soil depths after 100 years of paddy and non-paddy management. PR100: 100-year-old paddy soil; NPR100: 100-year-old non-paddy soil.
Scientific RepoRts | 7: 10818 | DOI:10.1038/s41598-017-10071-0 Changes in phosphorus speciation with prolonged management. The tidal wetland (TW) sediment, i.e. the parent material of soil formation for paddy and non-paddy soils, had a total P concentration of 760 mg P kg −1 in surface soil (Table 3). After land embankment, soil formation started with the salt marsh (SM) having a comparable P concentration (730 mg P kg −1 ) in the surface horizon. With prolonged arable management, total P concentrations increased, reaching on average significantly greater concentrations in the non-paddy topsoils (average of 1128 ± 198 mg P kg −1 ) than in those under paddy management (average of 890 ± 149 mg P kg −1 ; Table 3; p < 0.001). The total P recoveries in NaOH-Na 2 EDTA extracts of the youngest soil were extremely low (4.2-5.2% of total P in TW and SM respectively, Table 3). With prolonged management, the E P concentrations from paddy topsoil increased during pedogenesis to 12.0-31.4% of total P (Table 3). This gain in both absolute concentrations and relative proportions of E P occurred for both P i and P o species, although with different curve progression (Table 3, Fig. 1). Notably, P o concentrations increased only during early stages of soil formation and then remained relatively constant in older soils, both under paddy and non-paddy management (Table 3, Fig. 1b). Kinetic modeling revealed that accumulation of P o in paddy topsoils needed about 50 years more time to reach steady-state (Xe = 178 years) than in non-paddy topsoils, but differences in rate constants (k) were not significant when including the uncertainty of k into the assessment of Xe (Table 4 and Fig. S1).
Multiplying the 31 P-NMR results (in percent) by those for E P and factoring in changes in bulk density gave results in kg P ha −1 (Table 5), which showed that concentrations of extracted orthophosphate in paddy soil significantly increased (p < 0.001) with increasing duration of paddy rice cultivation. The concentrations of all other P compounds followed a similar trend through early stage of pedogenesis, i.e. being low or undetectable at the onset of soil development in TW and SM, but increasing thereafter in the paddy soils rapidly to a maximum that was obtained after no more than 300 years (Table 5). At the latter stage of the chronosequence, the concentrations of some P compounds like pyrophosphate, polyphosphate, orthophosphate monoesters, myo-IHP, other orthophosphate monoesters (calculated as total P in the orthophosphate monoester region after subtracting myo-IHP, scyllo-IHP, αand β-glycerophosphate and the mononucleotides) or glycerophosphates + mononucleotides even tended to decline, whereas those of other compounds could be sustained even after prolonged land use (e.g., scyllo-IHP, orthophosphate diesters, phosphonates; Table 5). In the non-paddy topsoils, the concentrations of all P compounds increased rapidly after 50 years of management, and then continued to increase slowly (except for orthophosphate and DNA) with subsequent management (up to 700 years; Table 5).
To better understand changes in P dynamics, we normalized the respective concentrations to the total amount of P i or P o extracted, i.e., we evaluated to which degree a given P i or P o compound contributed to overall P i or P o from 31 P NMR spectra. The results showed that paddy topsoils showed increasing proportions of total Pi for orthophosphate and declining ones for pyrophosphate during soil development, whereas non-paddy topsoils showed the opposite tendency ( Fig. S1a and b). The contributions of polyphosphates to extractable P i did not exhibit a clear temporal trend in both land-use systems, though proportions were larger in the paddy soils (Fig. S1c). In contrast to the P i forms, changes in P o composition were more consistent: the proportions of orthophosphate monoesters to extracted P o increased until a steady-state maximum was reached after 113 years of paddy and 91 years of non-paddy management ( Table 4, Fig. 2a). Moreover, orthophosphate monoesters constituted the highest proportion of total P o in both paddy and non-paddy topsoils at steady-state, though with higher proportions in the non-paddy soils (Fig. 2a). The myo-IHP and scyllo-IHP as typical orthophosphate monoesters showed similar accumulation patterns to total orthophosphate monoester, but they constituted only a small proportion of total P o in both paddy and non-paddy topsoils ( Fig. 2b and c, Table 5). The orthophosphate monoesters, myo-IHP, and scyllo-IHP were found to need a longer time to reach steady-state in paddy topsoils than in non-paddy topsoils (Table 4).
Although orthophosphate diesters may have a different origin than orthophosphate monoesters, the contribution of orthophosphate diesters to total P o also increased rapidly at early stages of pedogenesis until a steady-state plateau was reached in the paddy soils; no such trend was observed for the non-paddy topsoils (Fig. 2d). The sum of the diester degradation products (α/β-glycerophosphates from phospholipids and mononucleotides from RNA) also showed this curve progression to a plateau with higher proportions in the paddy topsoils (Fig. 2e). The DNA content was larger in the paddy topsoils than the non-paddy ones, and tended to increase with the increasing duration of management and related pedogenesis in the paddy soils (Fig. 2f). Only the phosphonates did not confirm a P accrual with time -their proportions declined with time in both paddy and non-paddy soils (Fig. 2g). To determine if the modelling results were affected by the longer time scale for the paddy soils, we redid our analyses using 700 years as the oldest date for both soil types. The trends observed when the 1000 and 2000  Table 3. Bulk density and concentrations of total P (P-total) and NaOH-Na 2 EDTA extracted P (E P ), inorganic (E p -P i ) and organic P (E p -P o ) concentration (kg P ha −1 and mg P kg −1 ) as well as their respective contributions to P-total or to E-P (%) in soils during 2000 years paddy and 700 years non-paddy management (n ± SD). TW: tidal wetland; SM: salt marsh; PR50-2000: 50-2000 years paddy soil; NPR50-700: 50-700 years non-paddy soil. years of non-paddy managements determined in NaOH-Na 2 EDTA extracts by solution 31 P-NMR spectroscopy (n ± SE). Mono-exponential regression (P < 0.05); no significant curvature is indicated by dashed line. Note: at t = 0 year the location was a tidal wetland and after 30 years it was still a salt marsh.
year old sites were included for the paddy soils were also observed with the analyzes up to 700 years only ( Fig. S2 and Table S2). The overall changes in organic P compounds with time are shown in Fig. 3.

Discussion
Changes in soil properties with time. We did not analyze any physical properties for these samples, nor any chemical properties beyond P forms by P-NMR, because this information is available in other studies 13,15,28 .
We have summarized the information here that is expected to be relevant for understanding P cycling in these soils. The pH values were 8.2 and 7.8 in the youngest soil (TW and SM; Table 6). With prolonged soil development, the pH values declined to 5.1 after 2000 years of paddy-rice management. This increasing acidification was related to the loss of carbonates (Table 6) 14 . In this regard, long-term paddy cultivation with submergence and drainage, as well as ploughing and puddling, resulted in faster decalcification than in non-paddy soils 14 , and thus in accelerated loss of related primary materials (e.g. CaCO 3 ), clays, and P-sorbing Fe oxides 3 . These processes likely led to the increase in P in non-paddy soils compared with soils under paddy management (Table 3), despite greater SOM accrual in the latter 13,14 . Higher yields of paddy rice than non-paddy crops also contributed to lower P storage in paddy soils. Concentrations of calcium (Ca), magnesium (Mg), Fe, manganese (Mn), and aluminum (Al) decreased with cultivation time in the paddy topsoils, so that the greatest concentrations of these elements were found mainly in the younger surface soils 15 . The dithionite-citrate-bicarbonate-(DCB-) extractable Fe in paddy and non-paddy chronosequences did not show large differences. However, greater concentrations of oxalate-extractable Fe (i.e. poorly crystalline Fe) in the paddy topsoils than in the non-paddy soils supported previous findings of enhanced redox cycling and accelerated weathering ( Table 6). The greater proportion of poorly crystalline Fe oxides was probably also responsible for the large proportion of mineral-associated soil organic matter in paddy soils 29 . The N and microbial residue accumulation were limited to initial stages of paddy soil development and restricted to the surface horizons 13 . Consequently, we assumed similar impacts on related P forms.
Phosphorus extraction and identification of phosphorus forms. Overall, NaOH-Na 2 EDTA extracted 12-31% of P-total in the paddy and non-paddy soils. We tested extraction methods to try to improve P recovery for these soils (see supplemental information); however, we found that a sequential repeated extraction did not increase the proportion of extracted P o and caused hydrolysis for some orthophosphate diesters (data not shown). The low P extraction yield in these soils may be related to the marine sediment parent materials and inputs of riverine sediments; most labile P was probably removed by sea water, leaving the majority of P occluded within Ca minerals. Low P recoveries (16-29% of total P) have been reported for NaOH-EDTA extracts of Ca-rich sea sediments 30, 31 , and similar low recoveries of P from Ca-rich soils have been reported for NaOH-Na 2 EDTA extracts 32, 33

Parameters
Unit k ± SE X 0 X e  Table 4. Kinetic parameters of the mono-exponential model (Eqn. 1: and for NaOH-NaF extracts 34 . The total P recoveries in NaOH-Na 2 EDTA extracts of the youngest soils without pronounced pedogenesis were extremely low (4.2-5.2% of total P in TW and SM respectively, Table 3). Turner et al. 22 suggested that low proportions of extracted P indicate that the majority of phosphate is present in primary minerals. The low TP and P o concentrations of the tidal wetland and salt marsh samples are consistent with other reports of P concentrations in the sediments of this region 35 . It is important to note, however, that inputs of riverine sediments to coastal areas in this region have been significantly altered by upstream water diversion over time 36 . As such, current nutrient concentrations in the tidal wetland and salt marsh may not fully reflect concentrations 700 or 2000 years ago, but should reflect conditions for the more recent soils in these chronosequence.

Organic phosphorus storage in paddy and non-paddy soils.
The paddy and non-paddy soils of the chronosequence were formed from marine sediment with an extremely low P o concentration, as indicated by the TW E P -P o concentrations (Table 3). This suggests that the majority of P o in the older paddy and non-paddy soils did not originate from marine P. After land embankment, soil formation led to an initial accumulation of P o from the few grasses and shrub vegetation growing in the topsoils 37 . With the start of cultivation, the P o stocks continued to increase exponentially (Fig. 1b), probably due to increased biomass inputs from the crops and manure applications 13 and from decreasing pH and changes in mineralogy 17 . The decomposition of organic matter is less efficient under anaerobic conditions during anthropogenic submergence of paddy topsoils. This results in a larger accumulation of organic matter relative to the non-paddy topsoils 14,27 . The gains in organic matter and the use of organic fertilizers such as manures likely promoted P o accrual in the surface soils 15 . Dense plough pans, developed from anthropogenic ploughing and puddling in the paddy soils 13 , reduced infiltration rates and thus blocked the vertical transport of organic matter into paddy subsoil 38,39 . As a result, the concentrations of all P o species generally declined more strongly with depth in the paddy soil profile than in the adjacent non-paddy counterpart (calculated by using Tables 1 and 2). We believe that this tendency is consistent with soil profiles of other ages, given that the 100-year-old paddy topsoil had some of the largest P o concentrations among the soils samples along the chronosequence (Table 3). Thus, the increase in paddy P o stocks over time with cultivation seems to be mainly driven by the accumulation of P in the topsoil but not by processes storing P in the subsoils. In addition, advanced Fe reduction in the paddy subsoil 15 likely hindered the effective binding of P in that part of the soil. The accumulation of P o in both land-use systems followed exponential growth curves to a maximum, indicating a saturation of P o concentrations and stocks in both systems over time. Such saturation implies both limited sorption capacity and steady-state conditions, i.e., at given amounts of P o stored, P o input and output rates must be similar. The findings are consistent with patterns of organic N accumulation in topsoils of both land-use systems, which reached plateaus within a similar time scales [78-193 years for N stock 13    Inorganic and organic phosphorus pools. Phosphorus speciation along a chronosequence is mainly regulated by changes in soil mineralogy for P stabilization, P inputs and biological production and utilization of various P compounds 2 . Alkali-extractable P i was classified as phosphate that is mainly associated with Fe and Al in soil 40,41 . It is a dynamic P pool and may act as an important P source for rice growth 42 . As P i concentrations continuously increased with increasing duration of paddy management (Table 5), our results implied a sustained change in soil nutrient composition for at least 2000 years. The results are different from native ecosystems where fertilization is lacking: results from a coastal dune chronosequence and nearby Franz Josef glacial chronosequence in New Zealand indicated that the concentrations of alkali-extractable phosphate increased only in the early stages of pedogenesis and then declined due to depletion of total P 2, 5 . The concentrations of pyrophosphate and polyphosphate were so small that they hardly contributed to the overall contents of total P i . The proportions of pyrophosphate even decreased with time in paddy soil, suggesting that it is not a stable P form in these environments. Low pyrophosphate proportions were also found in humic acid extracts of paddy soils in Philippines 43 . We conclude that pyrophosphate is easily decomposed, even when microbial activity may be reduced during the wet season under submerged anaerobic conditions in paddy soil 43 . Certainly, these changes in redox conditions promoted greater steady-state percentages of orthophosphate diesters in paddy than in the non-paddy soils (Fig. 2d). Orthophosphate diesters may accumulate due to reduced decomposition or increased microbial synthesis. The exact source assignment of orthophosphate diesters to plants and microbes is frequently ambiguous. However, paddy soils also contained larger contents of microbial N residues than their non-paddy counterparts 13 ; Therefore, it seems reasonable to assume co-accumulation of some microbial P forms that were rich in orthophosphate diesters 5 .
Intriguingly, and unlike other P o species, the portions of DNA continued to increase in paddy soils with time (Fig. 2f). The DNA extracted for 31 P-NMR may be present in the soil in microorganisms or plant residues, or may be stabilized on soil minerals, and stabilization in soil or in soil organic matter may be favored at lower soil pH 2,17 . Similarly, proportions of DNA increased with time in the 120,000 year Franz Josef post-glacial chronosequence, New Zealand, correlating significantly and positively with soil total organic carbon and N and negatively with soil pH 2 . Thus, with respect to accumulation of DNA with time, both the anthropogenically-modified ecosystems of this study and the native ecosystems in the New Zealand showed similar stabilization patterns of DNA with time.
The proportion of the sum of α/β-glycerophosphate and mononucleotides to total P o increased in the first 50 years and then slightly declined with prolonged cropping (Fig. 2e), which caused a similar tendency for total orthophosphate diesters. However, the tendency is equivocal because the highest proportion at 50 years is only sustained by one data point and may thus also be an outlier in our chronosequence; additionally, no other P species showed a similar tendency. Therefore, we ignored the highest point and considered the α/β-glycerophosphate and mononucleotides to have reached a plateau after steady-state conditions. Soil phospholipids and RNA are associated with living soil microbial biomass 2,[44][45][46] . In this line, a higher proportion of α/β-glycerophosphate and mononucleotides is consistent with increased microbial residues in the paddy soils compared with non-paddy soils that were observed in previous studies from these sites 13,15 .
Unlike the forest/shrub soil chronosequences in New Zealand with abundant inositol hexakisphosphate 1, 2, 5 , myo-IHP and scyllo-IHP constituted only a small proportion of orthophosphate monoesters in our study. In general, inositol phosphates may be stabilized by amorphous Fe oxides 2 , and close relationships between the contents of inositol phosphates and amorphous metal oxides have been found in a range of ecosystems 2, 47-50 . However, the submergence during paddy cultivation induced anaerobic conditions, under which amorphous Fe oxides are easily dissolved. Hence, associated inositol phosphates may be released or even lost with subsequent drainage or surface runoff. Additionally, there could be potentially larger inputs of inositol hexakisphosphate from plant materials (e.g. seeds) into the forest/shrub soil chronosequences in New Zealand than into the studied paddy soils 2 .

Site
Depth (cm) Horizon pH* Ca* (mg g −1 ) IC* a (mg g −1 ) Fe DCB * (mg g −1 ) Fe ox * (mg g −1 ) TW(n = 1) 2-30 8. Kinetic response of organic phosphorus to prolonged paddy and non-paddy management. The Walker & Syers model 51 of P transformation predicts that with increasing time of pedogenesis there will be a loss in total P, a depletion of primary mineral P pools, but an accumulation of P o stocks in early stage of soil development followed by a slow, subsequent decline with time. Our sites are different to those studies by Walker and Syers 51 in that our soils are fertilized to supplement P losses; besides, our land-use sequence covers a much shorter time scale to the chronosequence of Walker and Syers 51 . Our detailed assessment of P o species using 31 P NMR, including individual P o compounds and functional groups along 2000 years of rice paddy and non-paddy soil development demonstrated that: (i) the chemical nature of P o changes during pedogenesis, and (ii) soil paddy management strongly influences P o composition. All P o species accumulated rapidly in the early stages of pedogenesis. Thereafter, a plateau that is typical for steady-state conditions was reached. Studies in natural ecosystems 2, 6, 7 noted that P o stocks that accumulated during early soil development stages would decline gradually in later stages of the ecosystem, as a response to nutrient limitation. For example, soil P o accumulated under N limitation but declined under P limitation in natural ecosystems 2 . In our arable study sites, concentrations of available N were low in the young paddy systems 13 , so the initial gain in P o may be explained by similar processes as described by Turner et al. 2 ; in addition, the mere accrual of organic matter 14 co-accumulated P o . In older systems, however, nutrient limitations have been counterbalanced by fertilization, so P o losses in response to P limitation do not occur because these 'old' systems do not become P limited.

Materials and Methods
Samples. The  , as well as adjacent replicated sites under nonpaddy cultivation for 50, 100, 300, and 700 years (NPR50-700). Additionally, the soils with 100 years paddy and non-paddy managements (i.e. PR100 and NPR100) were sampled by soil horizon to around 100 cm for elucidating changes in P o composition with soil depth; for the other soils we only analyzed the top A horizons (~0-15 cm depth). More details about these soil samples are shown in Table 6 and may also be obtained from Roth et al. 13 and Lehndorff et al. 28 . It should be noted that the traditional paddy management practice in this region in China is a crop rotation of rice in the wet season followed by wheat or other upland crops in the dry season 53 . Total soil P concentrations were determined by melting a mixture of soils and lithium borate at 1000 °C for 30 min, followed by inductively coupled plasma optical emission spectrometry (ICP-OES) 54 .
Organic P analyses. Solution 31 P-NMR spectroscopy of alkaline extracts is commonly considered to be a reliable method for the quantification of P o 40 . Soil samples from each of the replicated plots at each location were extracted by shaking 4 g of air-dried soil with 40 mL of a solution containing 0.25 M NaOH and 0.05 M Na 2 EDTA for 4 h, followed by centrifugation at 10,000 × g for 30 min 23,55 . A 2-mL aliquot of each supernatant was used to determine total extracted P by ICP-mass spectrometer (MS); the remaining solutions were frozen at −80 °C, lyophilized, and ground. Each freeze-dried extract (~100 mg) was re-dissolved in 0.1 mL of deuterium oxide and 0.9 mL of a solution containing 1.0 M NaOH and 0.1 M Na 2 EDTA, and then immediately transferred to a 5-mm NMR tube. Solution 31 P-NMR spectra were obtained using a Bruker 600-MHz spectrometer equipped with a prodigy-probe (a broadband CryoProbe that uses N-cooled RF coils and preamplifiers to deliver a sensitivity enhancement over room temperature probes by a factor of 2 to 3 for nuclei from 15 N to 31 P). Extracts were measured with a D 2 O-field lock at room temperature, and chemical shifts were referenced to 85% orthophosphoric acid (0 ppm). The NMR parameters generally used were: 32 K data points, 2.6 s repetition delay, 0.7 s acquisition time, 30° pulse width and 10,000 scans. The delay time used here allows sufficient spin-lattice relaxation between scans for P compounds in NaOH-Na 2 EDTA (see supplemental information Fig. S3). Peak areas were calculated by integration on spectra processed with 7 and 2 Hz line-broadening, using NUTS software (2000 edition; Acorn NMR, Livermore, CA) and manual calculation. Individual P compounds were identified by manual inspection and the peak-picking subroutine in the NUTS software, based on their chemical shifts from reports in the literature 56 and by spiking selected samples with myo-inositol hexakisphosphate (myo-IHP) 1 . The orthophosphate peak in each spectrum was standardized to 6.0 ppm during processing 57,58 . Identified compounds and compound classes include orthophosphate (6 ppm), pyrophosphate (~−5 ppm), polyphosphate (−4 to −5, −5 to −50 ppm), orthophosphate monoesters (3 to 6, 6 to 7 ppm), orthophosphate diesters (3 to −4 ppm), and phosphonates (7 to 50 ppm). Because αand β-glycerophosphates and mononucleotides, which result from degradation of orthophosphate diesters during 31 P-NMR analysis 59 , were detected in orthophosphate monoester region, we assigned these compounds to orthophosphate diesters rather than to monoesters 58,60 . The concentrations of individual P species were calculated by multiplying 31 P-NMR proportions by the total NaOH-Na 2 EDTA extractable P concentration. The percentage of each P compound or compound class was calculated manually by integration across the entire spectrum and also by integrating smaller regions within each spectrum. The high-molecular weight (HMW) P o recently identified by McLaren et al. 61 is included with the unknown peaks in the P-mono-other category (Table S3). In our experience (B. Cade-Menun, unpublished data), peaks for these unknown HMW compounds are better resolved when lyophilized samples are prepared for NMR with NaOH-EDTA than Scientific RepoRts | 7: 10818 | DOI:10.1038/s41598-017-10071-0 with only water and D 2 O 61 , which is also reflected in the clear separation of all the myo-IHP peaks from the orthophosphate peak (Fig. S4).
Statistical evaluation. The total P concentrations between paddy and non-paddy soils were tested for significant differences (set to p < 0.05) by t test, one-way ANOVA was used to test significant differences of P compound concentrations among soil chronosequences, and regression functions of concentrations and distributions of various P i and P o during pedogenesis were calculated. All statistical analyses were conducted using Sigmaplot 12.5 for Windows. To estimate accumulation rates, a mono exponential regression model was used: where X t is the parameter of issue at cultivation time t (years), X e is the parameter stock at absolute equilibrium, X 0 is the initial parameter concentration in the tidal wetland (t = 0), and k is a rate constant 13 .
Data availability. All data generated or analysed during this study are included in this published article (and its Supplementary Information files).