Continuous decrease in soil organic matter despite increased plant productivity in an 80-years-old phosphorus-addition experiment

The study addresses how eight decades of agricultural management, including tillage and phosphorus (P) fertilization, affected soil organic matter (SOM) stocks, stoichiometry, P pools, and carbon (C) dynamics in a temperate cropland. The context is the need for sustainable agricultural practices that maintain soil fertility; however, most prior work assessed either soil element cycling or yields over relatively short periods. Long-term experiments are rare but critical to understand cumulative effects of tillage, fertilization, and biomass removal. Previous observations from long-term cropland trials showed mixed outcomes: some reported declines in soil total organic carbon (TOC) regardless of fertilization due to land-use change and aggregate disruption by plowing, while others found TOC increases with inorganic nutrient additions due to greater plant inputs. The impact of inorganic P fertilization on TOC is ambivalent: it can boost productivity and inputs, but P can also enhance SOM decomposition by desorbing organics from minerals. The objective was to quantify long-term changes in TOC, total nitrogen (TN), total organic P (TOP), soil P pools, and yields across control and two P addition treatments, and to evaluate SOM turnover using Δ14C measurements and modeling to infer transit times and pool dynamics.

The paper reviews long-term experiments showing divergent TOC trajectories under different managements. In Switzerland (Zurich Organic Fertilization Experiment) and Woburn (UK), TOC decreased across treatments, attributed to conversion to cropland and repeated plowing that disrupts aggregates and accelerates decomposition. Conversely, at Bad Lauchstädt (Germany) TOC increased under inorganic nutrient additions but decreased without nutrients, likely due to differences in plant-derived organic matter inputs. P fertilization often increases productivity and may raise TOC, particularly when combined with nitrogen (e.g., Bad Lauchstädt; La Estanzuela, Uruguay). However, some trials showed no TOC response to P fertilization (Zurich), and others reported that inorganic P can increase SOM decomposition or deplete TOC under N limitation, potentially via phosphate-induced desorption of adsorbed organic matter. Long-term P additions can alter plant-available P and organic P pools, including phytate, though stabilization mechanisms and soil mineralogy (e.g., high clay vs. Fe/Al oxides) modulate outcomes. This mixed evidence motivated the current long-term assessment integrating yields, SOM stocks, P pools, and radiocarbon-based C dynamics.

Study site and design: The long-term field experiment is at Lanna field station (58°20′49.9"N, 13°07′36.1"E), South-West Sweden (75 m a.s.l.), on Quaternary glacial clay overlain by ~30 cm postglacial silty clay loam. Mean annual temperature/precipitation were 6.1 °C/558 mm (1961–1990) and 7.3 °C/584 mm (1991–2020). Soil: Udertic Haploboroll; upper 0–20 cm is fine-textured with ~42–45% clay and ~15% fine silt; average pH 6.3 over 53 years. The experiment has two halves (started 1936 and 1941), with each treatment replicated twice per half (n=4 total). Treatments analyzed: Control (no P), P (1.75 g P m−2 yr−1 as superphosphate annually), and P6 (10.5 g P m−2 every 6th year; same average annual rate). All plots receive Ca(NO3)2 annually. Plot size is 7 m × 7 m. Management included an initial 7-year rotation with manure applications (2.0 kg m−2 applied twice per rotation in the first ~20 years); since late 1950s rotations shifted to mainly cereals, no manure thereafter, and soils are plowed annually to 20 cm with residues retained. Data collection: Yields were measured almost annually since inception. Topsoil (0–20 cm) samples were collected every six years before P application; samples were air-dried, sieved (<2 mm), and roots removed. Archived composite samples (per treatment) exist for 1968, 1977, 1980, 1989, 2001, 2003, 2013, 2015, 2019 (1968 only Control and P). In 2021, new samples were collected from all four replicates per treatment: five cores per plot (0–20 cm) composited within plot, air-dried, sieved, and analyzed separately. The time series thus spans 1968–2021 with 10 points. Grain P content was determined annually for 1995–2009. Chemical and isotope analyses: Plant-available P (P-AL) was extracted using acid ammonium lactate (0.1 M NH4-lactate + 0.4 M acetic acid, pH 3.75) with 1.5 h shaking, filtered (0.2 μm), and measured by ICP-OES (Avio 200, Perkin Elmer). Soil pH: 6.0 g soil in 18 mL H2O. Total organic P (TOP) was determined by ignition method (difference between inorganic P in 0.5 M H2SO4 extracts of ignited vs. non-ignited aliquots) with phosphate quantified by molybdenum blue (Murphy-Riley) on a continuous flow system (AA500, Seal). Total P (TP) was taken as P in the non-ignited extract. Phytate-P was isolated after extraction in 0.25 M NaOH + 50 mM EDTA followed by hypobromite oxidation; P in purified phytate solution measured by ICP-OES. Total organic C (TOC) and total N (TN) were measured by combustion (LECO CNS-2000). Inorganic C was also checked to confirm carbonate-free soils. Radiocarbon (Δ14C) of TOC was measured by MICADAS AMS (Ionplus) at the Max Planck Institute for Biogeochemistry, reported as Δ14C (‰) normalized for δ13C and decay-corrected. Other measurements: Particle size distribution was determined after H2O2 pretreatment by wet sieving and sedimentation. Grain P content was measured by boiling milled grains in 65% HNO3 and analyzing digest P by ICP-OES. Calculations and statistics: Stocks for 0–20 cm were calculated assuming a bulk density of 1.3 g cm−3 (constant over time). Relative yield of P and P6 treatments was computed yearly as (treatment yield / control yield) × 100. Treatment differences in 2021 were tested by ANOVA with Tukey’s post hoc, using Shapiro–Wilk and Levene’s tests to verify normality and homoscedasticity; significance threshold P < 0.05. Temporal trends were modeled using linear and monoexponential models; model fits compared via AIC and residual standard error. Analyses were conducted in R 4.1.1. Modeling: A two-pool C model (fast pool Cf and slow pool Cs) with series connection was used to simulate TOC stocks and Δ14C dynamics, driven by annual yields and the atmospheric bomb 14C curve. Equations: dCf/dt = γY − kf Cf − αf Cf; dCs/dt = αf Cf − ks Cs; with initial conditions Cf(0) = β C0, Cs(0) = (1 − β) C0, where C0 is initial TOC. Two input scenarios were tested: (1) constant C input equal to 150% of long-term average yield C (per Bolinder et al.), and (2) variable C input proportional to annual yields with proportion γ. Parameters (γ, kf, ks, αf, β) were estimated via inverse modeling using SoilR for radiocarbon-enabled simulations and FME for optimization: initial classical optimization to define priors followed by Markov chain Monte Carlo to obtain posterior distributions. The best-fit parameters and uncertainties were used to compute mean and median transit time distributions following Metzler & Sierra.

• TOC, TN, and TOP in the upper 0–20 cm declined linearly over five decades by approximately 13.7–13.8% for TOC and TN and 11.6% for TOP (1968–2021), irrespective of P fertilization. Mean annual TOC loss across treatments was 0.14 t C ha−1 yr−1. • Crop yields increased from about 220 g dry weight m−2 in 1936 to more than 500 g m−2 in the 2010s. Phosphorus addition increased yields by 26–30% (P and P6 vs. Control) over 1936–2021. • Soil TP and TOP stocks in 0–20 cm did not differ significantly among Control and P treatments after ~80 years, indicating that plant P uptake in the Control was largely from subsoil (>20 cm). P-AL stocks were significantly higher in P-fertilized treatments (P < 0.01), while phytate-P and TOP did not differ significantly by treatment in 2021. • Grain P concentration averaged 3.7 mg g−1; average yield (1936–2021) was 385 g m−2 yr−1, implying mean P removal of ~1.42 g P m−2 yr−1 by grain harvest. In the last decade, yields >600 g m−2 yr−1 (especially with P fertilization) removed >2.2 g P m−2 yr−1, exceeding the annual P addition rate (1.75 g P m−2 yr−1), implying subsoil P supply. • Δ14C of TOC was very low in 1968 (≈ −120‰), remained below 0‰ through 2021, yet showed evidence of bomb 14C incorporation; this pattern reflects large stocks of older C and relatively low capacity for Δ14C enrichment in this cold, clay-rich soil. • Two-pool model fits reproduced TOC and Δ14C dynamics under both constant and yield-proportional input scenarios. Estimated γ (variable-input model) was 1.443 ± 0.48, indicating C inputs ≈144% of crop yield C. Despite differences in pool sizes and turnover between scenarios, both required an overall increase in decomposition rate to reconcile TOC decline with bomb 14C incorporation. • Mean transit time of C was below 10 years (8.72 years for constant-input model; 4.86 years for variable-input model), and median transit times were below 7.5 years (5.86 and 3.29 years, respectively), indicating that a large fraction of new C inputs are respired rapidly (about 50% within 2–6 years), with only small proportions persisting beyond a decade. • Decades of P addition did not significantly alter TOC stocks, TOP, Δ14C, or SOM C:P stoichiometry, suggesting limited impact of P fertilization on SOM stabilization in this clay-rich soil where desorption-driven SOM loss appears minor.

The findings show that, despite substantial increases in plant productivity over eight decades and significant yield gains with P fertilization, SOM in the topsoil has continuously declined. This indicates that the current management—annual plowing to 20 cm, residue retention without organic amendments, and mineral fertilization—enhanced decomposition relative to inputs, likely through aggregate disruption and accelerated mineralization. Modeling of TOC and Δ14C constrained C cycling, indicating short mean and median transit times and rapid respiration of new inputs. The necessity to invoke increased decomposition rates in both modeling scenarios points to a long-term shift in SOM turnover dynamics rather than solely changes in inputs. The absence of treatment differences in TP and TOP in the top 20 cm, alongside elevated P-AL in fertilized plots, suggests that high yields in the Control were sustained by P uptake from subsoil reserves (>20 cm), challenging assumptions that topsoil P dominates plant nutrition in such systems. The lack of TOC response to P fertilization and stable SOM C:P ratios are consistent with strong sorptive stabilization of phosphorylated organic compounds by clay minerals, mitigating P-induced SOM desorption effects observed in Fe/Al oxide-rich soils. Overall, the results address the central question by demonstrating that long-term tillage and management can deplete SOM stocks even under increased productivity, highlighting a sustainability gap: soil fertility and other properties (e.g., aggregation, water retention) may be at risk if SOM decline continues.

High crop yields in the control treatment were supported by plant P uptake from below 20 cm, as evidenced by similar topsoil TP and TOP stocks across treatments after ~80 years. TOC and TOP stocks in the top 20 cm decreased by 13.8% and 11.6%, respectively, over 53 years across all treatments, indicating that current agricultural practice at this long-term site is not sustainable due to continuous SOM depletion. Radiocarbon-informed modeling showed mean C transit times below 10 years and rapid respiration of new inputs, implying that TOC responds slowly to changes in organic C inputs. The study underscores the importance of long-term (>50 years) agricultural experiments for understanding plant–soil element cycling and for guiding sustainable management that maintains SOM.

Stock calculations assumed a constant bulk density (1.3 g cm−3), not accounting for potential temporal changes due to management. Soil sampling and analyses focused on the uppermost 20 cm; subsoil pools and dynamics, though implicated in P uptake, were not directly quantified. Archived samples prior to 2019 were composited by treatment, limiting within-treatment replication for those years. The Δ14C and TOC time series comprise 10 time points from 1968–2021, with missing data for some treatments in 1968, which may constrain temporal resolution. Modeling results showed uncertainty in pool partitioning and turnover (different pool decline patterns under constant vs. variable input scenarios), reflecting parameter uncertainty and structural assumptions about inputs and transfers.