Subsoils—a sink for excess fertilizer P but a minor contribution to P plant nutrition: evidence from long-term fertilization trials

The phosphorus (P) stocks of arable subsoils not only influence crop production but also fertilizer P sequestration. However, the extent of this influence is largely unknown. This study aimed to (i) determine the extent of P sequestration with soil depth, (ii) analyze P speciation after long-term P fertilization, and (iii) compare soil P tests in predicting crop yields. We analyzed four long-term fertilizer trials in Germany to a depth of 90 cm. Treatments received either mineral or organic P, or a combination of both, for 16 to 113 years. We determined inorganic and organic P pools using sequential extraction, and P speciation using 31P nuclear magnetic resonance (NMR) and X-ray absorption near edge structure (XANES) spectroscopy. In addition, we applied three P soil tests, double-lactate (DL), calcium acetate lactate (CAL), and diffusive gradients in thin films (DGT). The results suggested that plants are capable of mobilizing P from deeper soil layers when there is a negative P budget of the topsoil. However, fertilization mostly only showed insignificant effects on P pools, which were most pronounced in the topsoil, with a 1.6- to 4.4-fold increase in labile inorganic P (Pi; resin-P, NaHCO3–Pi) after mineral fertilization and a 0- to 1.9-fold increase of organic P (Po; NaHCO3–Po, NaOH–Po) after organic P fertilization. The differences in Po and Pi speciation were mainly controlled by site-specific factors, e.g., soil properties or soil management practice rather than by fertilization. When modeling crop yield response using the Mitscherlich equation, we obtained the highest R2 (R2 = 0.61, P < 0.001) among the soil P tests when using topsoil PDGT. However, the fit became less pronounced when incorporating the subsoil. We conclude that if the soil has a good P supply, the majority of P taken up by plants originates from the topsoil and that the DGT method is a mechanistic surrogate of P plant uptake. Thus, DGT is a basis for optimization of P fertilizer recommendation to add as much P fertilizer as required to sustain crop yields but as low as necessary to prevent harmful P leaching of excess fertilizer P.


Background
Phosphorus (P) is indispensable for plant growth and one of the essential nutrient elements for agriculture. Despite being present in soils in large quantities, not only in topsoils but also in the subsoil, most P is not readily available chemically to plants [33,49], often making P fertilization essential to maintaining high crop yields. However, there are uncertainties regarding finite P rock reserves [26], discussions about an uneven distribution of deposits [8], and environmental concerns such as the eutrophication of water bodies due to surplus P application [61]. Siebers et al. Environ Sci Eur (2021) 33:60 Thus, there is a need for a more sufficient (not more than required) and efficient (higher yields per applied P) use of P in the future. This includes, for example, the recycling of P from waste and residues, as well as increasing utilization of existing P in the soil through optimal fertilization and plant management.
The subsoil (here defined as soil underneath the plough layer, approximately > 30 cm) can serve as a P reservoir for plant nutrition, as it can store substantial amounts of P ranging from 25 to 70% of total P in the soil profile [6,7,29,32]. Indeed, 10% to 50% of subsoil P is acquired by crops when topsoil P is low [35]. However, to be readily available for plant nutrition, the P must be in a suitable chemical form and physically accessible to plants. Koch et al. [32] found that subsoil P species in an agriculturally used Stagnic Cambisol were dominated by P associated with iron (Fe) and aluminum (Al), irrespective of fertilizer treatment. In a subsequent study, they estimated how much of these subsoil P species contribute to plant nutrition and found that up to 30% of the plant P requirements were covered by these P species under optimal soil physical conditions in the early plant growth phase [31]. The subsoil not only provides P to the plant, it also serves as a sink for excess fertilizer P above plant requirements. The degree of P saturation (DPS) is an indicator for the risk of P losses. Soils with DPS values ranging from 25 to 40% have a high risk of P loss, either by surface runoff or by leaching [48]. Thus, to optimize P fertilization and crop management, detailed knowledge about availability and accessibility is essential not only for topsoil but also for subsoil P resources.
In particular, the use of agricultural waste and residues for P fertilization, such as compost, manure, or green manure, is of importance for an efficient and sustainable closed-loop production system. Long-term field experiments are most suitable for examining P fertilization management on the soil P status. They enable us to study complex fertilizer turnover processes in soil operating on different time scales and thus provide an overview of the effectiveness of fertilizer management on nutrient mobilization, transformation, translocation, and uptake [22,30,76,79].
To obtain a detailed overview of the soil P status and the dynamics of topsoil and subsoils, the application of complementary methods is essential. A combination of wet-chemical extractions and spectroscopic approaches to determine P pools and their speciation has seen widespread application in agriculturally used soils under different fertilization management [9, 22,32,34,70]. A promising method to determine labile P is the diffusive gradients in thin films (DGT) technique [18,33]. The DGT approach mimics the actions of the plant roots in accessing available P at the soil-solution interface. This is achieved through a local decrease in the P concentration of the soil solution via sorption to an infinite iron-based P sink. It thus simultaneously triggers the resupply of P by diffusion from the bulk solution and desorption from the solid soil phase. A 33 P isotope dilution study showed that among eight agricultural soil P tests, P concentrations estimated by DGT exhibited the closest relation to the P uptake by plant roots [59]. This test is independent of the soil pH or calcium carbonate content [40], which may influence the erratic results of other extractionbased tests, such as double-lactate (DL) or calcium acetate lactate (CAL) [77]. Besides, if P is the only limiting nutrient element and if P uptake is diffusion-driven, DGT appears to accurately determine yield response to P fertilization [41,47,62]. Therefore, we assume that the DGT approach is more accurate in estimating crop yields than liquid chemical extraction-based procedures.
The DGT approach, in combination with sequential P fractionation methods [27], provides a detailed overview of P pools with different availabilities. However, as these methods are operationally defined, no precise information is obtained about the chemical P speciation in the soil. This can be overcome using spectroscopic approaches. Commonly applied spectroscopic approaches here are P K-edge X-ray absorption near edge structure (XANES) spectroscopy as well as 31 P liquid nuclear magnetic resonance ( 31 P NMR) spectroscopy. P XANES spectroscopy is an element-specific and non-invasive direct P speciation method that is mainly sensitive to the speciation of inorganic P (P i ) species [9, 32,33]. However, 31 P NMR enables the identification and speciation of organic P (P o ) forms [11,12,33,67]. Thus, a combination of these wet chemical and spectroscopic techniques provides a comprehensive overview of the availability and chemical form of P in soils.
There are numerous studies on the effect of different P fertilizers on P forms in the soil [9,20,32,36,38]. However, to the best of our knowledge, none of these studies comprises a broad spectrum of P fertilizer treatments together with analyses of the subsoil of different soil types. Here, we included three long-term field experiments that exhibit different soils, analyzed the topsoil and subsoils, and applied a multi-method approach including wet-chemical and spectroscopic approaches. This helps to provide an all-encompassing picture about how the effects of fertilization, soil type, and management on P pools are also relevant for crop yield prediction. To the best of our knowledge, this approach has never been applied to such an extent. With this study, we aim (i) to assess P availability and speciation in soils after long-term application at different soil depths and (ii) to evaluate the suitability of the DGT approach in predicting crop yields compared to other soil P tests. We studied four long-term field experiments with different durations of P fertilizer management, ranging from 16 up to 113 years up to the time of sampling at three sites. All experimental sites had the same P fertilizer treatments, namely no P addition, mineral P addition, organic P addition in the form of compost, green manure or farmyard manure, and a combined mineral and organic treatment. We sampled all treatments down to a depth of 90 cm. Topsoil and subsoil samples were subjected to a combination of DGT P analysis, sequential fractionation, P K-edge XANES spectroscopy, and 31 P NMR spectroscopy to demonstrate the long-term effects of fertilizer P management on P pools and species.

Experimental sites and soil sampling
We collected soil samples from four different agricultural long-term experimental sites for basic characteristics, sequential P fractionation, P DGT analyses, P XANES spectroscopy, and 31 P NMR spectroscopy. However, the agricultural long-term experimental site at the University of Rostock was analyzed for P DGT only, since data on sequential P fractionation, 31 P XANES, and NMR spectroscopy were previously published by Koch et al. [32]. Details on every experimental site can be found in Additional file 1: Table S1. The Dahlem and Bad Lauchstädt sites both received farmyard manure with substantial inorganic P contents [22]. However, to make a distinction from P applied by mineral fertilizers, the farmyard manure treatments are abbreviated as "organic P" or "mineral + organic P" (for a combined farmyard manure and mineral fertilizer treatment) throughout the manuscript.
Three soil cores were taken with a mechanical soil auger (diameter 3 cm) randomly distributed over each plot down to a depth of 90 cm. The soil cores obtained from Dahlem were divided into four segments, representing depths of 0-30 cm, 30-50 cm, 50-70 cm, and 70-90 cm, respectively, and slightly adjusted to horizons if required. The soil profiles were as follows: Ap (0-30 cm depth), EBw (30-50 cm depth), BtE (50-70 cm depth), and Bt (70-90 cm depth). The soil cores from Dürnast Freising and Bad Lauchstädt were divided into three segments, representing depths of 0-30 cm, 30-60 cm, and 60-90 cm, respectively. As no field replicates were available in Bad Lauchstädt, three soil cores from one plot were analyzed separately. The Chernozem soil at the Bad Lauchstädt sampling site has a very thick A horizon of about 50-60 cm [3]. The soil profile was as follows: former Ap (0-30 cm depth), Ah (30-45 cm), Ck + Ahk (45-60 cm), Ck (60-90 cm). We combined the Ahk and Ck + Ahk horizons and defined them as an unplowed A horizon. To enable direct comparisons among all three field experiments, we define the upper plowed 0-30 cm layer as topsoil and the unplowed 30-60 cm layer as well as the deepest layer (60-90 cm) as subsoil. For Dahlem and Freising Dürnast, soil core segments from each plot were well-mixed on-site and sub-sampled for the present study. All samples were air-dried and ground to pass through a 2-mm sieve.

General soil analyses
The soil texture was determined according to DIN ISO 11,277 (< 63 µm to 2 µm silt, < 2 µm clay). Soil pH (0.01 M CaCl 2 ) was measured according to Altermann et al. [3]. Soil carbon (C) and nitrogen (N) analyses were performed using the dry combustion method followed by thermal conductivity detection of the released trace gases (vario MICRO cube, Elementar, Hanau, Germany). Total elemental concentrations of P from bulk soil (P tot ), iron (Fe tot ), aluminum (Al tot ), manganese (Mn tot ), magnesium (Mg tot ), potassium (K tot ), calcium (Ca tot ), and sodium (Na tot ) were determined after the microwaveassisted digestion of 150 mg with 0.7 mL HNO 3 and 2 mL HCl using an inductively coupled plasma-optical emission spectrometer (ICP-OES; Thermo Fisher iCAP ™ 7600). The bulk densities of the samples were calculated for each depth using the soil dry weights after drying at 105 °C and the auger volume at specific sampling depths. Total elemental stocks were calculated for each sampling depth with a respective bulk density (Additional file 1: Table S1). In addition, ammonium oxalate extractable P (P OX ), Fe (Fe OX ), and Al (Al OX ) were determined according to Schwertmann [58], DL extractable P was determined according to Riehm [53,71], and CAL extractable P was determined according to Schüller [57,72]. The total P concentrations in the extracts were determined by ICP-OES. Using the P OX , Fe OX , and Al OX , the DPS was calculated as described in Koch et al. [32]. (Additional file 1: Fig. S1; see Additional file 1 for DPS results). All chemicals used for general soil analyses at least had the purity level for analysis.

Soil analyses for P pools
To test the desorptive P release from soils, we used the DGT approach described in the literature [18,78]. The DGT sampler (DGT Research Ltd, Lancaster, UK) consisted of a hydrogel layer containing an ion resin that is covered by a diffusion layer and protected by a membrane. It is a solute sampling technique, which can be applied to bulk soil estimations of the plant-available P [41,59,65]. We deployed the DGT device to the soil brought to 100% water holding capacity for 24 h at 21 °C. After deployment, the DGT units were dismantled and the ferrihydrite gel eluted in 1 mL of 1 M HCl solution for at least 24 h. The concentration of total P in the eluent was measured by ICP-OES. Standard DGT formulae described by Zhang et al. [78] were used to calculate the P DGT concentrations.
All soil samples were also subjected to sequential P fractionation as described in [32]. In total, 0.5 g of the soil sample was sequentially extracted using (1) H 2 O + resin, (2) 0.5 M NaHCO 3 , (3) 0.1 M NaOH, and (4) 1 M H 2 SO 4 following the decreasing extractability of P. The molybdate-reactive P [hereafter referred to as inorganic P (P i )] concentration of the filtrates was determined photometrically using the molybdenum blue method [46]. Total P was determined by ICP-OES and molybdate unreactive P, principally organically complexed P [63] [hereafter referred to as organic P (P o )] was calculated as the difference between total P and inorganic P. All chemicals used for soil analyses of P pools at least had the purity level for analysis.

Soil analysis for P speciation
Due to the relatively low P concentrations, especially in the subsoils that were not sufficient to yield P K-edge XANES spectra with sufficient signal-to-noise ratio, all samples were relatively enriched in P by isolating the < 20 µm size fraction (in agreement with [32], since most P is concentrated in the clay-to medium-silt-sized fraction [37]. In this way, we achieved a P tot enrichment by a factor of up to 6 compared with the bulk soil samples. This ensured a sufficient signal-to-noise ratio for all samples. P K-edge XANES spectra were collected using the Soft X-ray Microcharacterization Beamline (SXRMB) at the Canadian Light Source (CLS) in Saskatoon, Canada. Due to limited beamtime, only samples from the control and combined mineral + organic fertilization treatments were analyzed using P K-edge XANES spectroscopy. The air-dried and homogenized samples were spread as a thin film onto a double-sided carbon tape and mounted onto a copper sample holder before being placed in the vacuum chamber. The spectra were recorded in fluorescence yield mode (samples) and total electron yield mode (reference standards), respectively, at photon energies between 2121 and 2215 eV. For each sample, at least two scans were recorded. Subsequent data treatment and evaluation, such as averaging of raw spectra (improving signal-tonoise ratio), background correction (pre-edge range: − 20 to 6 eV relative to E0), and normalization (normalization range: + 28 eV to + 50 eV normalization order: 2), as well as linear combination fitting (LCF), were performed using the ATHENA software package (Demeter 0.9.26) [51]. Linear combination fitting was performed on averaged, normalized spectra in the energy range between − 12 and + 35 eV, relative to the E0. No floating of the E0 was allowed during fitting. Fitting was performed for all possible binary to quaternary combinations of the 11 most probable P reference standards in the sample: CaHPO 4 , CaHPO 4 × 2H 2 O, Ca(H 2 PO4) 2 × H2O, amorphous calcium phosphate (ACP), AlPO 4 , FePO 4 × 4H 2 O, P adsorbed on gibbsite, P adsorbed on goethite, and (C 6 H 18 O 24 P 6 ·Na·yH 2 O) = phytic acid sodium salt hydrate). The r-factor values were used as goodness-offit criteria. Significance between fits was evaluated using the Hamilton test (P < 0.05) [14] with the number of independent data points calculated by Athena, estimated as data range divided by core-hole lifetime broadening. If, despite the addition of a further reference compound to a fit, the r-factor was not significantly better according to the Hamilton test, the fit with a lower number of reference compounds was chosen. If fits with the same number of reference compounds were not significantly different from each other, obtained proportions of the reference compounds and the r-factors of all significantly different fits were averaged. All chemicals used for XANES spectroscopy at least had the purity level for analysis. An exception to this was the standard P adsorbed on gibbsite and P adsorbed on goethite, as these compounds were produced by others [2,25].
For 31 P NMR spectroscopy, the same treatments as for P K-edge XANES spectroscopy were analyzed. Sample preparation for the 31 P NMR measurement followed a procedure previously described in Cade-Menun et al. [13] and Turner [66]: 2 g of bulk soil were shaken for 4 h in 40 mL of 0.25 M NaOH-0.05 M EDTA mixture at room temperature. After centrifugation at 14,000×g for 30 min, the supernatant was frozen and lyophilized. P recovery after extraction was summarized in Table 1. In total, 100 mg of freeze-dried solids was redissolved in 500 µL of a mixture of NaOD and D 2 O with a pH of 13. For quantification, 100 µL of methylenediphosphonic acid solution (MDPA in NaOD/D 2 O, 0.8 mg mL −1 ) was added as an internal reference standard [10]. The sample was centrifuged again at 14,000×g for 30 min and the supernatant was decanted into a 5-mm NMR tube. The NMR spectra were acquired on a Bruker 600 MHz spectrometer equipped with a CryoProbe. Acquisition parameters were listed as follows: 0.55 acquisition time, 13,700 scans, inverse-gated proton decoupling, and 293.15 K. Considering the sample was excited by a 30° pulse calibrated at 4 μs, a 3-time T1 delay was required between each scan to ensure complete relaxation of all P nuclei [10]. The T1 of the individual sample was estimated based on the P/(Fe + Mn) ratio [42]. For signal assignment, reference standards from Sigma Aldrich, including α-,β-glycerophosphate, and myo-inositol hexakisphosphate (myo-IHP) were spiked into soil extracts. To identify RNA hydrolysis products, an RNA reference standard was deliberately hydrolyzed in a NaOD/ D 2 O mixture and spiked into the soil extract. The NMR spectra were processed using MestReNova software (version 8.1.2-11880) with 2 Hz line broadening. All chemicals used for 31 P NMR spectroscopy at least had the purity level for analysis.

Calculation of P budgets
The P budgets were calculated as follows: The P input results from organic fertilization (P OF ) and mineral fertilization (P MF ), and P output results from P uptake by the harvested crop products (P Plant ). For the Dahlem experiment, the P budgets were calculated for 4 years from 2014 to 2017 using annual measurement data (crop yields, amount of farmyard manure, dry matter contents) and laboratory data (P contents of harvested products and farmyard manure). As no information was available regarding P stocks in the topsoil and subsoils from the beginning of the long-term experiments, we had to calculate the theoretical initial P stocks by adding the P budget (positive or negative) of the control treatment to measured P stocks of the control treatment. Using this value, we were then able to calculate the theoretical P stocks that should be present in the soil, taking into account the P budgets for the entire duration of the longterm experimental trial by offsetting the P stocks of the treatments estimated for different profile depths with the P budgets (Additional file 1: Table S2). We were then able to estimate the deviation of the theoretical P stocks from the measured P stocks, while incorporating different soil depths (Additional file 1: Table S2).

Yield response assessment
To evaluate the relationship between the recorded crop yield data of all sites and the available soil P as determined by different methods (P DL , P CAL , and P DGT ), we determined the maximum yield according to Mason et al. [41] by fitting the Mitscherlich equation: where y zero is the yield of the control and y zero + a is the maximum yield reached with the highest P application calculated using the model. The crop response to fertilization was determined by calculating the relative yield (%) using the maximum yield calculated with Eq. (1): We included all sites and treatments for crop response modeling using the Mitscherlich equation (Eq. 2) within (1) P budgets = P OF + P MF − P Plant .
(2) y = y zero + a 1 − e −bx , one fit [41]. To estimate the depth-dependent quality of the fit, we determined the Mitscherlich fit for three depth intervals separately, namely 0-30 cm, 0-60 cm, and 0-90 cm. In doing so, it is possible to estimate the contribution of the subsoil layers to crop prediction. The values of P DGT , P DL , and P CAL for Dahlem at depths of 30-50 cm and 50-70 cm were combined and included in the fits for the 0-60 cm depth of the other sites. We plotted soil P data (P DGT , P CAL , or P DL ) against relative yields from Eq. (3) and used the obtained curve to estimate fitting parameters y zero + a and b. By inserting these fitting parameters and soil P data (P DGT , P CAL , P DL ) into Eq.
(2), we then obtained relative yield data as estimated by the Mitscherlich fit, y fit . The quality of the fit was estimated by correlating the y fit values with the values obtained through Eq. (3). Using this approach, it is then possible to evaluate the performance of different soil P tests to estimate crop response to fertilization under field conditions.

Statistical analyses
Statistical analyses were performed using IBM SPSS 22 software (IBM SPSS Statistics V. 22.0, 2013, IBM Inc.). We tested total elemental concentrations, sequentially extracted P fractions, as well as P CAL , P DL , and P DGT for normal distribution using the Shapiro-Wilk test (P < 0.05) and for homogeneity of variances using the Brown-Forsythe test (P < 0.05). We considered samples from different depths as dependent samples. To account for this, we performed a 2-way repeated measures ANOVA, with fertilizer treatment as the 1st factor and depth intervals-treated as repeated measures-as the 2nd factor. If significant differences occurred, we used the Tukey HSD test for post hoc separation of means (P < 0.05). In addition, simple linear regression was performed correlating the values of P DL or P CAL with P DGT . We determined correlations separately for the depth steps 0-30 cm, 0-60/70 cm, and 0-90 cm analogous to the approach of the Mitscherlich fit described above. No statistical analyses were conducted for samples from Bad Lauchstädt, as no field replicates were available; therefore, differences between treatments and depths were pointed out as trends.

P budgets
When comparing the measured P stocks with theoretical P stocks, the deviation was greatest in Bad Lauchstädt when only incorporating the topsoil (0-30 cm) with an underestimation of P stocks ranging from 12 to 46% (Additional file 1: Table S2). This deviation became lower by incorporating the other depth into the calculations and the smallest deviation was estimated by incorporating the complete soil profile from 0 to 90 cm (Additional file 1: Table S2). However, in Rostock, no downward movement of P was observed [32] and thus the deviation was lowest when only considering the first 30 cm of the soil profile (Additional file 1: Table S2). In Dahlem, there was no difference in deviation when incorporating the subsoil layers, except for the organic P treatment.

P fractions
Total P stocks mostly did not reflect P fertilization in Dahlem and Freising Dürnast, while in Bad Lauchstädt, total P stock in the topsoil tended to increase after mineral P and combined mineral + organic fertilization compared to the control (Additional file 1: Table S3). In comparison with the control, organic fertilization tended to increase not only the total P stocks of the topsoil but also the total P stocks in subsoils in Bad Lauchstädt. To highlight P pools of different availabilities in relation to total P stocks, we performed a sequential extraction.
Hedley extractability (sum of extracted P during sequential Hedley fractionation in relation to total P) ranged from 18 to 53% for Dahlem, 60% to 101% for Dürnast, and 64% to 90% for Bad Lauchstädt. In general, stocks of stable P pools (H 2 SO 4 -P and residual-P) were greatest among all pools for all sites and treatments, followed by labile (resin-P, NaHCO 3 -P) and moderately labile (NaOH-P) P pools (Additional file 1: Table S4, Fig. 1).
With increasing soil depth, the labile and moderately labile stocks decreased, which to some extent was also reflected in their proportions on total P stocks (Additional file 1: Table S4, Fig. 1). The proportions of stable P (H 2 SO 4 -P, Dahlem; H 2 SO 4 -P and residual-P (organic and mineral + organic P treatments), Bad Lauchstädt) tended to increase in deeper soil layers compared to the topsoil, but was only significant in a few cases (Additional file 1: Table S4). The treatment effects were less pronounced than the depth effects. There was a significant increase of NaHCO 3 -P i stocks after mineral P and combined mineral + organic P fertilization in the topsoil compared to the control, while NaOH-P o stocks significantly increased in the topsoil after organic P fertilization in Dahlem. This increase in the P o of the topsoil after organic fertilization was also visible in Freising Dürnast for NaOH-P o , although the increase was not significant. These trends were visible in stocks as well as in proportions. In Bad Lauchstädt, treatments mostly affected labile inorganic P i (resin-P, NaHCO 3 -P i , NaOH-P i ), raising their stocks as well as proportions compared to the control for all depth layers (Additional file 1: Table S4, Fig. 1). NaOH-P o (stocks and proportions) in Bad Lauchstädt also tended to increase in subsoil with combined mineral + organic fertilization. In addition, organic P fertilization tended to increase the H 2 SO 4 -P concentrations and proportions of the topsoil in Bad Lauchstädt. In summary, with the labile (resin-P, NaHCO 3 -P), moderately labile (NaOH-P), and stable (H 2 SO 4 -P and residual-P) P pools, similar effects were visible as for single fraction stocks and proportions (Fig. 1).
In particular, the impacts on the labile P pool, which is known to reflect plant-available P to a certain extent, is of interest after P fertilization. For this reason, we also determined P DGT , which was highest in the topsoil for all sites, irrespective of fertilization treatment (Fig. 2). In addition, between treatments, significant effects were most prominent in the topsoil and the upper subsoil layer, except for Freising Dürnast, which only showed significant effects in the topsoil. The control treatment receiving no P fertilization was mostly lowest in P DGT in Dahlem for the first and second soil depths, while there was no significant difference between the control and the fertilization treatments for the other sites. However, although not statistically tested, Bad Lauchstädt also exhibited the lowest P DGT in the control treatment for the first and second soil depths (Fig. 2). In addition, in Dahlem, the highest P DGT concentration was found for the mineral P treatment in the topsoil (0-30 cm) and the upper subsoil layer (30-50 cm), while for the other subsoil layers, no significant differences were observed. The only significant treatment effect on P DGT in Freising Dürnast was that topsoil P DGT was the lowest among all treatments by a significant margin for organic P fertilization, including the control without P fertilization. For Bad Lauchstädt, no statistical analyses were possible. However, it seems that similar to Dahlem, the control treatment exhibited the lowest P DGT concentration and the mineral + organic P treatment the highest P DGT concentration in the two topsoil layers (0-30 cm and 30-60 cm). Similar findings were also made for the Rostock site with significantly higher P DGT concentrations in the mineral + organic P treatment in the topsoil (0-30 cm) and upper subsoil layer (30-60 cm), while all other treatments were not significantly different.

Solution 31 P NMR spectroscopy
The soil P o speciation was revealed by 31 P NMR analyses after extraction with the NaOH-EDTA mixture. P recovery after extraction, which appears to be higher in topsoils than subsoils, ranged between 10 and 64% relative to P t stock (Table 1). In the monoester region, the signals of α-,β-glycerophosphate, and myo-IHP were identified by spiking experiments, while the spiking of RNA hydrolysis products did not correspond to any signals (Additional file 1: Figs. S2, S3). The signals at 4.1 ppm and 6.8 ppm were assigned to scyllo-IHP and isomerized IHP, respectively [68,69]. Freising Dürnast and Dahlem soils showed similar monoester composition, with all predominant Fig. 1 Stock of labile (Resin-P, NaHCO 3 -P i,o ), moderately labile (NaOH-P i,o ), and stable (H 2 SO 4 -P t and Residual-P t ) P fractions. Samples were obtained from the long-term experimental sites Dahlem, Freising Dürnast, and Bad Lauchstädt at different soil depths. Error bars represent standard deviations of the mean (n = 3) signals from myo-and scyllo-IHP. In Bad Lauchstädt, the myo-IHP signals were weak but a strong unknown signal was observed at 4.4 ppm. The signals of α-,βglycerophosphate were low and considered as hydrolyzed phospholipids for quantification (Table 1).
A trend of decreasing P recovery with increasing sampling depth was observed, especially for Dahlem and Bad Lauchstädt. The soil P pool at the Dahlem site was mostly inorganic, with orthophosphate accounting for 80-100% of P t . For the other two sites, the proportion of orthophosphate, although lower than Dahlem, was still the most abundant P form (39-82%) followed by monoesters (17-59%). Phosphonate, diesters, and pyrophosphate were only detected in small proportions and were not present in all samples ( In all samples, the contents of orthophosphate and extractable P o were highest in the topsoil and gradually decreased with depth (with a few exceptions at Freising Dürnast). The content of pyrophosphate was also higher in the topsoil than in the layers below. The proportions of orthophosphate and extractable P o only showed a clear depth-dependent trend at Dahlem; no such trend was observed at Freising Dürnast and Bad Lauchstädt.
Compared to the control plots, the extractable P t content in both topsoils and subsoils increased after the application of fertilizer, especially in Bad Lauchstädt soils where the most pronounced P enrichment effect was observed. In addition, the growth in extractable P t stocks after fertilizer treatment was mostly attributed to orthophosphate, since only the orthophosphate proportion increased, while the share of other P categories such as monoesters declined. At Bad Lauchstädt, fertilization decreased the P o content in the topsoil, but the opposite effect was observed in subsoils. At the other two sites, the effect of fertilization on total P o content was unclear.

P-XANES analyses
The LCF of the P K-edge spectra (Additional file 1: Fig.  S10) indicated that the Dahlem and Freising Dürnast sites were dominated by Fe-associated P (26-73%), while soils in Bad Lauchstädt were dominated by Ca-associated P (27-56%). The proportions of organic-associated P were lowest in Dahlem (3-10%), whereas significant proportions were found for Freising Dürnast (15-27%) and Bad Lauchstädt (17-27%) ( Table 2). It should be noted that focusing on the < 20 µm size fraction might also have resulted in a shift in the P i /P o ratio due to a higher relative enrichment of P o compounds, which are known to associate and accumulate in the clay-to the medium-sized silt fraction [37]. Treatment effects were not consistent among sites. No clear treatment effect was Fig. 2 Plant available P concentration as determined by the diffusive gradients in thin films (DGT) approach. Samples were obtained from the long-term experimental sites Dahlem, Freising Dürnast, Bad Lauchstädt, and Rostock at different soil depths. Different lowercase letters indicate significant differences among treatments within each soil depth (P < 0.05). Error bars represent standard deviations of the mean (n = 3). No statistical analysis could be conducted for samples from Bad Lauchstädt as no field replicates were available observed for Freising Dürnast. By contrast, in Dahlem, mineral + organic P fertilization resulted in a relative decrease in Fe-associated P in favor of Al in the topsoil, whereas organic-associated P remained the same. In the subsoil, P speciation was very similar to the control and a possible treatment effect was only indicated by an increase in organic-associated P in both subsoil layers. A decrease in Fe-associated P in favor of Al-associated P relative to the control was also observed in the topsoil in Bad Lauchstädt. However, in contrast to Dahlem, organic-associated P here decreased in the topsoil whereas Ca-associated P increased. In the subsoil, mineral and organic P fertilization only affected proportions in Bad Lauchstädt, with generally decreasing proportions of P associated with Ca-P and increasing proportions of organic-associated P at a depth of 60-90 cm, while it remained similar at a depth of 30-60 cm (Table 2).

Crop yield response
The stocks of the DL extractable P (P DL ) and CAL extractable P (P CAL ) (Additional file 1: Table S5) for all three sites (Dahlem, Freising Dürnast, Bad Lauchstädt) followed the same trend as P DGT (Additional file 1: Table S5, Fig. 2). When correlating P DL or P CAL with P DGT significant positive linear relationships were found for both P DL and P CAL for the individual depths of 0-30 cm, 0-60/70 cm, and 0-90 cm with R 2 ranging between 0.27 and 0.60. The coefficient of determination was throughout higher for P CAL than for P DL (Fig. 3). Incorporating the subsoil had no effect on P DL and P DGT regression, but resulted in higher R 2 for P CAL and P DGT regressions.
Rostock, Dahlem, and Bad Lauchstädt were responsive to P fertilization, i.e. crop yields were significantly higher for the highest P fertilizer application level compared to the control treatment (Additional file 1: Table S6). Therefore, using Eq. 2, maximum yields were estimated at 119 dt ha −1 (maize), 7 dt ha −1 (winter rye), and 60 dt ha −1 (spring barley) for Rostock, Dahlem, and Bad Lauchstädt, respectively. Freising Dürnast did not show significantly higher crop yields with the highest P application level (Additional file 1: Table S6). We therefore determined the maximum yield through the application of all P levels except the control treatment. We thus obtained a mean value of 127 dt ha −1 (maize). In doing so, we included all sites and treatments for crop response modeling using the Mitscherlich equation (Eq. 2) [41]. Crop response to Table 1 Depth resolved phosphorus (P) stocks and the composition of NaOH-EDTA extractable soil P, as assessed by the liquid 31

P-NMR spectroscopy of homogenized composite samples (n = 3) of three long-term experimental sites
The values in parentheses (%) represent the individual P forms as a proportion of total NaOH-EDTA extractable soil P. The content of diesters and monoesters were corrected by considering P hydrolysis during extraction

Management
Depth P recovery Phosphonate Orthophophate Monoesters Diesters Pyrophosphate Extractable P o Extractable P t > 20 ppm 6.7 to 6.0 ppm 6.9 to 6.7, 6.0     applied fertilizer P was modeled well, including all sites and treatments, using the Mitscherlich relationship for three soil P tests for the topsoil, whereas fitting became worse when subsequently adding deeper soil layers to the fit (Fig. 4, Additional file 1: Table S7). The best and significant coefficients of determination for the regression fits R 2 were obtained for P DGT , followed by P CAL , while fits for P DL were the worst (Fig. 4).

Influence of fertilization on P stocks and pools
The increasing P stock with fertilization in Bad Lauchstädt was also found by others [44]. The increase was anticipated, as there was a positive P budget calculated for the mineral (1299 kg ha −1 ), organic (214 kg ha −1 ), and combined mineral + organic fertilization (2921 kg ha −1 ), whereas the control treatment had a negative P budget (− 1569 kg ha −1 ) until the time of soil sampling [44,75]. In contrast, in Freising Dürnast, unchanged topsoil P stocks after organic P fertilization were surprising, as the P budget for the control was − 1105 kg ha −1 and the budget of, for example, the mineral P treatment was 324 kg ha −1 for an experimental duration of 37 years. Significantly increased P stocks had been anticipated for that treatment (Additional file 1: Table S3). Zicker et al. [79] attributed the gap in P budgets and the topsoil P stocks to the potential of the plant to also acquire P from deeper soil layers and the downward movement of fertilizer P from upper soil layers, which is assumed to apply here.
An uptake from deeper soil layers and especially a downward movement of fertilizer P was also indicated in the data of Freising Dürnast and even more pronounced in the data of Bad Lauchstädt for this study. When considering the DPS-indicating the risk of P losses-the topsoil DPS in particular was > 30% for the Dürnast and Bad Lauchstädt sites for treatments receiving P (Additional file 1: Fig. S1), indicating a risk of applied P leaching into subsoil layers. At both sites, organic as well as combined mineral + organic P fertilization tended to increase total P stocks (although not significantly) in the subsoils compared to the unfertilized control, which was most pronounced for the combined fertilization (Additional file 1: Table S3). In the literature, the same trends for the application of organic fertilizers were found for other long-term experimental sites [22,39]. The authors attributed the downward leaching of P after organic fertilization to the co-leaching of P complexing compounds or compounds released from organic fertilizer that compete with the organic fertilizer-derived P for sorption sites [28,52]. This also explains our observation that the combined mineral + organic P fertilization in particular facilitates a downward movement of fertilizer P. Among the three tested sites, the enrichment of P in the subsoils was highest in Bad Lauchstädt, followed by Freising Dürnast, whereas no enrichment was observed at Dahlem, despite the fact that DPS was > 40% (Additional file 1: Fig.  S1). This suggests that the extent of the downward movement and enrichment of P in the subsoil appears to be site-specific.
Possible site-specific factors are, for example, the amount of applied P vs. plant P uptake, the duration of the experiment, and soil properties. For this study, it is very likely that besides duration-which particularly increases the contrast between control and P treatments-the variations in soil properties between the sites are the main driving factors. The Retisol of Dahlem with the clay-enriched argic subsoil horizon may have limited a downward movement of P into deeper subsoil layers. By contrast, in the Cambisol in Freising Dürnast, and especially the Chernozem in Bad Lauchstädt (well-structured, and pore-rich chernic horizon down to 60 cm), the downward movement of water is facilitated.
If these considerations are correct, not only the topsoil but also the subsoil must be considered when calculating the P budgets. The greatest deviation from theoretical P stocks found for Bad Lauchstädt when only incorporating the topsoil (0-30 cm) (Additional file 1: Table S2) is a result of the downward movement of P due to the chernic horizon down to 60 cm. In Dahlem, a downward movement of P was mostly prevented due to the clay-enriched subsoil horizons starting at a depth of ~ 60 cm. Therefore, no distinct trend for the extent of deviation was observed when incorporating different soil layers. Another important factor is the existing duration of the fertilization regime. Prior to the time of sampling, the experimental site Bad Lauchstädt had been managed for 113 years, whereas the site in Rostock was only established 16 years ago, explaining the lowest deviations when only considering the topsoil. The factors soil properties and time were also reflected when comparing Rostock with Freising Dürnast. Both soils are Cambisols but the deviation (See figure on next page.) Fig. 3 Relationship between P DL and P CAL concentrations and P DGT concentrations. Analyses were conducted for all treatments (control, mineral P, organic P, mineral + organic P), all experimental sites (Dahlem, Freising Dürnast, Bad Lauchstädt, Rostock), and the sampling depths 0-30 cm, 0-60/70 cm, and 0-90 cm. Significant correlations were highlighted with *(P < 0.05), **(P < 0.01), and ***(P < 0.001). DL double-lactate, CAL calcium acetate lactate, DGT diffusive gradients in thin films Siebers et al. Environ Sci Eur (2021) 33:60 in Freising Dürnast remained the same or even became 0 for the combined mineral + organic P treatment when incorporating the complete soil profile sampled, most likely an effect of the longer duration of this experiment (16 vs. 37 years). These results demonstrate that for P budget calculations, several parameters have to be taken into account when evaluating P fertilization strategies. Fertilization mostly increased the proportion of P stored in the topsoil. The P fertilization not only affected the distribution of total P stocks along the soil profile, but also different P pools to a certain extent. For example, the increase of NaOH-P o in the topsoil after organic P fertilization and partial increase of H 2 SO 4 -P (Bad Lauchstädt) (Additional file 1: Table S4) was anticipated due to the introduction of organic materials also incorporating occluded P i and P o , i.e. stable humic acids and larger humus complexes or monoesters that are covered or physically encapsulated [5,17,32,73]. However, a change in the P fractions of the subsoil was not very pronounced. Mineral P fertilizers (alone or in combination with organic fertilizers) affected P i fractions in all soil depths, but were always most pronounced in the topsoil. This was mainly visible in the labile P pool (H 2 O-P i + NaHCO 3 -P i ), but it was also reflected in P DGT at the Bad Lauchstädt and Rostock sites and in P DL and P CAL at Bad Lauchstädt for the mineral + organic treatment (Fig. 2, Additional file 1: Table S5). The fact that in Bad Lauchstädt, the mineral + organic P fertilizer was more effective in increasing labile P (labile P i , P DGT , P DL , and P CAL ) than the mineral P alone is likely due to the fertilization history of the plot with higher application levels until 1981 (Additional file 1: Table S1). This also applies to the Rostock site with P application levels above plant requirements for the combined mineral + organic P treatment [32].
Furthermore, a possible competition for sorption sites between fertilizer P i and organic anions and acids released during the decomposition of the organic fertilizers [32,52] increases the soluble soil P pool. There is also a clear difference between the types of organic fertilizer visible in P DGT concentrations. While the application of green manure or just compost (also containing substantial amounts of green waste; [32]) decreased the P DGT concentration of the topsoil by − 46% (Freising Dürnast) or only slightly increased it by 14% (Rostock), the application of farmyard manure substantially increased the P DGT concentrations of the topsoil by 30% in Dahlem and by 86% (0-60 cm) in Bad Lauchstädt. These results support findings made by other scientists, claiming that soils Fig. 4 Mitscherlich relationship of plant-available soil P (DL, CAL, and DGT) with crop response. n = 48; response was modeled using the Mitscherlich equation. Significant coefficients of determination (R 2 ) were highlighted with *(P < 0.05), **(P < 0.01), and ***(P < 0.001). DL double-lactate, CAL calcium acetate lactate, DGT diffusive gradients in thin films amended with manure exhibited higher plant-available P than soils receiving compost [22].
The sites receiving green manure or compost (Freising Dürnast and Rostock) both exhibited comparable P proportions distributed among P fractions, which can be attributed to both soils being Cambisols (Additional file 1: Table S4, [32]). For both sites, treatment effects were also less pronounced and mostly insignificant, including with respect to inorganic P fractions like P DGT . Koch et al. [32] concluded that 16 years of continuous fertilization in Rostock was not a sufficient period of time to lead to pronounced management-induced differences in P fractions, which might also apply to Freising Dürnast with 37 years of continuous fertilization management prior to the time of soil sampling.

Influence of fertilization on phosphorus speciation
Current P research still lacks a complete understanding on the extent to which the P forms present in the topsoil and subsoil are physically accessible and chemically available to the crop [29,35,64]. The use of soil depthresolved sequential P fractionation in combination with spectroscopic methods helped to show how P forms at the molecular level are influenced by fertilization management.
Using 31 P NMR spectroscopy, we found a decreasing P extractability (NaOH-EDTA extraction) with increasing soil depth (Table 1), which is consistent with previous studies [22,32,55]. The dominance of orthophosphate stocks, regardless of the soil type or the crop rotation, was in line with other studies [1,32]. Along the soil profile, the degree of the accumulation of NaOH-EDTA extractable P i and P o forms was site-specific. This can be explained by the differences in the amount of applied P and the differences between the duration of the field experiments, which both potentially intensify contrasts of the P fertilization effects in the soil profile [22]. For instance, relatively small changes in the NaOH-EDTA extractable P i and P o stocks along the soil profile relative to the control were found in Freising Dürnast, which corresponds to the shortest duration among all three sites until time of sampling (36 vs. 92 or 113 years).
In general, similar to labile P and P DGT , the orthophosphate stocks were the highest in topsoils and gradually decreased with depth ( Figs. 1 and 2, Table 1), thus confirming that there are larger stocks of labile P in the topsoils than in the subsoil. Furthermore, the increase in stocks of orthophosphate in the subsoil after fertilization treatment in relation to the control supports a downward movement of labile P i , which also applied for monoesters in the subsoil, especially at the Bad Lauchstädt site.
Monoesters in the topsoil at Bad Lauchstädt were hardly influenced by P fertilization, likely indicating a balance between the input, mineralization, and stabilization of NaOH-EDTA extractable monoesters after 113 years of management. Therefore, the retention capacity of topsoil for monoesters appears saturated, facilitating the leaching of surplus monoesters into the subsoil. However, our data suggest that the observed difference in the accumulation patterns of ortho-P monoesters along the soil profiles between Freising Dürnast, Bad Lauchstädt, and Dahlem was also influenced by the above-mentioned site-specific soil properties at Dahlem restricting the transport of P into deeper subsoil layers.
The abundance of diesters, even after including their hydrolysis products, was generally low in all soils. By contrast, monoesters (0-59%) were the dominant P o form over dieters (0-1.1%), which was in agreement with other studies on arable soils [16,32,43]. It reflected an enhanced mineralization of the labile diester-P due to cultivation practices. For example, [63] showed that cropping and fertilizer input limited P o diversity in temperate arable soils.
The difference in P o speciation in the topsoil was more pronounced between sites, rather than between treatments, suggesting that the type of P fertilization did not significantly affect the P o speciation in the soil. This is in line with previous studies on long-term fertilization trials [4,19,32,56] and implies that at all three sites, differences in P o speciation were mainly controlled by site-specific factors, such as soil properties or soil management practice. However, most of the above-mentioned studies were restricted to P o speciation of the topsoil. Results from this study indicated that the effect of treatment in the subsoil was not consistent over all sites, e.g., accumulation of monoesters in the Bad Lauchstädt subsoil. This confirmed that P in Bad Lauchstädt is more prone to downward movement because of the soil properties, but other factors, including trial duration, application levels, initial P sorption capacity, and DPS may also play a role here.
The proportion of P o detected by P XANES correlated significantly with the proportions of P o on total P extracted with NaOH-EDTA (R 2 = 0.73, P < 0.001) and extractable P o quantified by 31 P NMR (R 2 = 0.53, P < 0.001). This confirms a reasonable comparability of both methods despite the known limitation of P K-edge XANES in identifying P o in complex matrices like soils [33]. In agreement with sequential fractionation (H 2 SO 4 -P), Ca-P substantially contributed to P speciation at all three sites. In Bad Lauchstädt, P associated with Ca increased with depth up to 56% (Table 2). This was in line with the highest H 2 SO 4 -P stocks among all studied sites comprising of relatively stable P most likely associated with Ca and Mg minerals [15,23,50] and reflected the calcareous parent material of the Haplic Chernozem in Bad Lauchstädt [3]. Similarly, the majority of Ca-associated P being found in 60-90 cm subsoil of the control can be attributed to the presence of primary P minerals in the calcareous material in the deeper subsoil. The higher proportions of Ca-P in the unfertilized topsoil at Freising Dürnast and Dahlem compared to the underlying subsoil layer (30-60 cm) are most likely a result of the continuous lime application and accumulation of Ca-P in this otherwise rather decalcified part of the soil profile.
Treatment effects on P speciation were again inconsistent among sites. Combined mineral + organic P fertilization resulted in almost no difference in P speciation in Freising Dürnast, which was in agreement with 31 P NMR. This was likely due to the relatively short duration of the trial, with changes too small to be detected by P XANES spectroscopy. The increasing proportions of Al-P in the topsoils at the expense of Fe-P after the application of fertilizer in Dahlem and Bad Lauchstädt confirmed results from other scientists showing that mineral fertilizer P was first associated with Al compounds in the topsoil [20,21]. The fact that proportions of Al-P in the subsoils with or without fertilization showed no great difference suggests that the leached fertilizer P most likely accumulated in other P forms in the subsoil. The P forms were also determined by site-specific properties. For instance, at Dahlem, leached P accumulated in the upper subsoil layer in the form of Ca-P as a result of continuous liming and the high content of Ca-P in organic farmyard manures [22]. The accumulation of P o in the 30-50 cm subsoil layer was in agreement with the 31 P NMR results and can be attributed to an enhanced downward transport of fertilizer P due to the low P sorption capacities and high DPS (Additional file 1: Fig. S1) of the sandy and clay-depleted horizons above the argic horizon starting at soil depths of around 65 cm. By contrast, the observed increase in P o , especially in the 60-90 cm subsoil layer, suggests that the risk of P o leaching seams higher in Bad Lauchstädt compared to Dahlem. This is likely due to the less restricted downward P transport within the porerich chernic horizon and the higher past and current organic P application rates.

Relation of P DGT to other soil P tests and consequences for crop response prediction
When looking at the linear correlations of extracted soil P with P DGT , it is apparent that to some extent, all methods access a similar soil P pool commonly known as plant-available P [18,54]. Thus, the increase in labile P after fertilization, especially after the addition of mineral fertilizer, which was in agreement with other studies [7, 22,32], was not only reflected in sequentially extracted P pools, but also in P CAL , P DL , and P DGT (Additional file 1: Table S5; Fig. 2).
The P DGT concentration can be viewed as the highly labile fraction of P CAL or P DL (Fig. 3), which is replenished over time within the P CAL or P DL pool. Correlations were generally comparable or better for the P CAL than for P DL , which is likely a result of the soil pH. In general, CAL should be applied for soils exhibiting a pH > 6, while DL is used for soils with a pH < 6 to also include apatitic phosphate in acidic soils [57]. As the pH of most samples was > 6, CAL performed better in extracting the labile P pool and thus correlated better to P DGT . Higher coefficients of determination found for correlations between P CAL and P DGT when also incorporating the subsoils is most probably due to slightly higher pH values in subsoils compared to the topsoil. However, although correlations were highly significant, the R 2 is only between 0.27 and 0.60, suggesting that the results of these methods cannot be fully transferred to each other by simple regressions, as was the case in this study. Therefore, further research on a much larger and more diverse sample set is needed to test the usability of additional soil parameters in multiple regressions so as to further improve predictability and, more generally, the usability of such transfer functions.
We also determine the ability of the different soil P tests to model crop yield response using the Mitscherlich equation. The highest R 2 were obtained using the P DGT data, which is in accordance with literature claiming that DGT accurately determines yield response to P fertilization in pot experiments [74] and field trials [40,41,45,47,59,60,62]. Up till now, most studies that related P predictions to yield based on P DGT were performed on Australian or tropical topsoils and only one study examined agriculturally used soils in Europe [47]. However, within this extensive study, only the topsoils of eleven different experimental fields were considered and there was no focus on subsoils. To evaluate the subsoils for crop production in agriculturally used soil in Europe, crop prediction models should also include subsoil layers, with this study representing a first step to closing this gap.
The coefficients of determination for the regression fit R 2 predicting yield response to soil P of different soil P tests were consistently better in topsoil than when incorporating one or all subsoil layers (Fig. 4). This implies that the majority of P taken up by plants originates from the topsoil. However, since the inclusion of subsoil layers only resulted in a significant prediction of yield response for P DGT -even when gradually decreasing with depththis suggests that the plant at least partially absorbed P DGT from the subsoil. It is known that plans can also acquire P from the subsoil [7,29], but to what extent and under which possible plant-specific conditions remain unclear. It is believed that P deficiencies in the topsoil can stimulate P acquisition from the subsoil [29], but at the same time, only a sufficient overall nutrient status in the topsoil leads to a high plant P uptake, facilitating root growth and, in turn, the acquisition of P from the subsoil. For instance, [7] was able to show that under N-limiting conditions, the utilization of subsoil P resources was restricted.

Conclusion
In summary, the literature and our own results have shown that plants are also capable of mobilizing P from deeper soil layers when there is a negative P budget in the topsoil [24,32,79]. The results of this study represent a good basis for future calculations on the influence of different subsoil layers for P plant uptake, as we provide the initial values of the stocks down to 90 cm for four test sites. We conclude that not only the topsoil but also deeper soil layers, the duration of the experiment, and soil properties have to be considered to guarantee a balanced P budget and reduce the risks of P leaching. Although the leaching of fertilizer P down the complete soil profile (0-90 cm) was observed for some sites, the influence of different P fertilizers on P pools and speciation was not very strong and mostly management-and site-specific, as well as being most pronounced in the topsoil. The topsoil also mainly contributed to the plant P uptake, as reflected in the Mitscherlich fits.
Although subsoil stocks were larger (among all subsoil layers) and the available P stocks of the subsoil were comparable or slightly lower than the available P stocks of the topsoil, there was no need for the plant to acquire substantial amounts of subsoil P due to the good P supply of the topsoil. For example, even after almost 100 years without P fertilization in Dahlem, there was only around a 10% yield decrease of wheat in the control treatment, but no significant reduction in the labile P stocks of the subsoils. This also indicates that the accessibility for plants is lower, despite the subsoil volume that can potentially be rooted by plants being many times higher than the topsoil volume. It can therefore be deduced that P fertilization tends to fill up the P reservoirs of the subsoil-the extent to which is dependent on site properties and management history-rather than largely contributing to plant nutrition in soils with good P status. A generalization of fertilization effects over different sites is rather difficult, as they largely depend on soil properties and fertilization management. However, we were able to show that, depending on the type of organic P fertilizer, the proportion of labile P was affected differently, with farmyard manure P exhibiting higher plant availability than green manure or compost. Among the tested soil P tests, the DGT method was best at predicting the plant response to P application, regardless of soil properties and crop type. Thus, we conclude that the DGT method is a mechanistic surrogate of P plant uptake and superior to equilibrium-type chemical extraction procedures at estimating potential crop yields, as it is independent of factors like soil type or soil pH.
Additional file 1. Results on 'General soil parameters' and discussion on 'Influence of fertilization on C and N stocks' . Table S1. Details about long-term experimental fields. Table S2. Total P stocks and P budgets. Table S3. General soil characteristics. Table S4. Sequential P fractionation data. Table S5. P DL and P CAL stocks. Table S6. Crop yields. Table S7. Mitscherlich fitting parameters. Figure S1. DPS for long-term experimental sites. Figure S2 and Figure S3. NMR spectra of NMR spiking experiment. Figure S4 to Figure S9. NMR spectra of the soil samples. Figure  S10. P K-edge XANES spectra of the soil samples.