Soil phosphorus over a period of agricultural change in Scotland

Yield increases by use of phosphorus (P) fertiliser has been a defining aspect of last century's agriculture, but in many countries P inputs are now being regulated to avoid ecological damage and improve agricultural sustainability. We tested the hypothesis that P forms and functions in agricultural soils have changed over a period of 50–70 years of agricultural change using topsoils from 35 agricultural sites in North East Scotland compared at ‘original’ and ‘resampled’ timepoints. The original dried and archived samples were collected between 1951 and 1981, and resampling took place in the autumn of 2017 with all samples reanalysed at the same time and by same methods. A range of soil extractants, from labile (water soluble) P, carbon (C) and nitrogen (N) forms, to surface exchangeable P (acid ammonium oxalate extraction), to complexed P (NaOH‐EDTA extraction), were compared with agronomic soil test P (modified Morgan; MM‐P), soil pH, and, on a subset of five time pairs (n = 10), specific P compounds using 31P NMR spectroscopy. Whilst MM‐P had declined significantly in time, all soils were originally and remained at high soil P status. In contrast, no change in labile P pools, nor in the more refractory P pools that replenish labile P, were observed. In the subset analysed by NMR, orthophosphate P consistently increased in concentration. This suggests that soil P shows a slow (multidecadal) response to the documented reductions in fertiliser P for the study region. Significantly greater soil pH (by 0.4 units) over time and reduced water‐soluble C and N forms suggests that labile P in contemporary soils exists in a biogeochemical environment of enhanced solubility and in excess to other available macronutrients that control its microbial cycling. This study represents a long and important period of changing drivers acting on soil P change with potential to improve understanding future soil P trajectories. The value of archived historical soils samples is shown, if care in interpretation of storage implications is made.

changed over a period of 50-70 years of agricultural change using topsoils from 35 agricultural sites in North East Scotland compared at 'original' and 'resampled' timepoints. The original dried and archived samples were collected between 1951 and 1981, and resampling took place in the autumn of 2017 with all samples reanalysed at the same time and by same methods. A range of soil extractants, from labile (water soluble) P, carbon (C) and nitrogen (N) forms, to surface exchangeable P (acid ammonium oxalate extraction), to complexed P (NaOH-EDTA extraction), were compared with agronomic soil test P (modified Morgan; MM-P), soil pH, and, on a subset of five time pairs (n = 10), specific P compounds using 31 P NMR spectroscopy. Whilst MM-P had declined significantly in time, all soils were originally and remained at high soil P status. In contrast, no change in labile P pools, nor in the more refractory P pools that replenish labile P, were observed. In the subset analysed by NMR, orthophosphate P consistently increased in concentration. This suggests that soil P shows a slow (multidecadal) response to the documented reductions in fertiliser P for the study region. Significantly greater soil pH (by 0.4 units) over time and reduced water-soluble C and N forms suggests that labile P in contemporary soils exists in a biogeochemical environment of enhanced solubility and in excess to other available macronutrients that control its microbial cycling. This study represents a long and important period of changing drivers acting on soil P change with potential to improve understanding future soil P trajectories. The value of archived historical soils samples is shown, if care in interpretation of storage implications is made.

| INTRODUCTION
The use of phosphorus (P) fertiliser has been one of the defining contributors to productive agriculture since the Green Revolution during the middle of the last century (Elser & Haygarth, 2020;Jarvie et al., 2019;Yates & Boyd, 1965). However, these increased yields have come at the cost of dependency upon the declining resources of P rock reserves and eutrophication of water bodies downstream (George et al., 2016). In this context, it is important to understand the long-term effects of P fertiliser additions on agricultural soil nutrient chemistry from the mid-20th century to present day. Knowledge on how the present soil nutrient conditions came about is necessary to inform future soil status trajectories and improved soil management and agricultural practise (Johnston & Poulton, 2019).
Since the mid-20th century there have been complex technically-, policy-and societally-driven changes in UK agriculture as elsewhere around the world. Change in the United Kingdom, and North East (NE) Scotland region of the present study (summarised in Figure 1), have involved changes in land cover (grassland-arable, woodland-farmed land changes) and management (e.g. drainage, field size), F I G U R E 1 Summary timeline of grouped factors relevant to change in national (Scotland) and regional (NE Scotland) agricultural practises over the period of the sample timepoints. References are [1] National Statistics Scotland (2016); [2] Domburg et al. (1998); [3] Edwards et al. (2016); [4] Buckingham et al. (2014) cropping (e.g., use of winter vs. spring cereals), animal management, use of fertilisers and other agronomic inputs (e.g. liming). The last 40 years has seen only a small $1% area increase in Scotland's farmed area (National Statistics Scotland, 2016). In NE Scotland, areas of grassland, root fodder crops and oats declined since 1960s, with increasing areas of barley and adoption of oil seed rape and winter cereals (Domburg et al., 1998). Cattle numbers rose until the 1980s then declined in NE Scotland, balanced by increased numbers of sheep, pigs and poultry (Domburg et al., 1998), with continuing sharp declines in numbers of cattle locally in recent decades with change from grassland to arable land use (Dunn et al., 2014). Fertiliser usage, crop yield and nutrient offtake combine to affect net nutrient accumulation. A surplus of P has been widely demonstrated across UK agricultural soils, although this has been reported as also reducing since the 1970s . Fertiliser P usage in United Kingdom increased from the 1940s but has declined from start of the current century (Department for Environment, Food and Rural Affairs [Defra], 2018; Withers & Foy, 2001). Annual rates of inorganic P fertiliser application for UK arable and grassland were approximately 40 and 30 kg P/ha, respectively in the 1960s (Yates & Boyd, 1965), changing to 50 and 25 kg/ha in the 1980s, and dropping to 30 and 11 kg P/ha in 2018 (Defra, 2018). In Scotland, steady overall chemical P fertiliser usage from 1960 to 1990 had mixed trends across crop types. In NE Scotland, P inputs declined for grassland, root crops and cereals, except for steady high applications to potatoes and increased rates for winter cereals, whilst manure (mostly cattle) remained a steady part ($11 kg P/ha/year) of P inputs until the 1990s (Domburg et al., 1998). In contrast, nitrogen (N) fertiliser usage increased greatly in NE Scotland to all crops until the decade of decline from 2000 and now remains some 30% lower in Scotland than at peak (National Statistics Scotland, 2016). Fertiliser P reductions were driven by costs and regulation and advice in response to negative effects of P surplus, particularly eutrophication in downstream water courses from areas of P fertiliser usage by runoff and leaching (Haygarth et al., 2005). The soil P surplus has been termed 'legacy' P and while it represents ongoing environmental risks (Powers et al., 2016), especially in the face of increased diffuse pollution climate drivers (Ockenden et al., 2017), it also represents an opportunity as a store of potentially available P for agriculture (Haygarth et al., 2014).
Soil P management is crucial to balancing production and environmental goals. Soils comprise a myriad of chemical P forms, which may be grouped into pools relating to their function within the P cycle in soil, often defined by their availability to different extractants and analysis by various methods (Kruse et al., 2015). In agronomy, plant available P (often referred to as soil test P, e.g. Khan et al., 2018) is used to advise on fertiliser recommendations for soil crop combinations that aim to balance crop P availability against the risk of loss of P by water movement and subsequent negative environmental issues (Roberts & Johnston, 2015). Although an important part of effective agricultural practise, agronomic tests of soil P are not enough to enable a full picture of the cycling of P in the soil, even as it relates to agronomic purposes. Weak extractions like water can inform on mobile P available for leaching or biological processes, whereas stronger extractants increasingly yield more refractory P forms that are more complexed or strongly bound to surfaces, transformations of which serve to replenish immediately available fractions (Haygarth et al., 2013). Reversible transfer of inorganic P (P i ) between pools of varying plant availability is important (Johnston & Poulton, 2019). Organic P (P o ) is also known to have an important role in P cycling in soil, 31 P NMR being the primary method used to investigate the detail of P o forms (Cade-Menun, 2017).
Land use and management change can have a profound effect on the quantity of soil P and distribution of soil P between various functional fractions and extractable pools. Arable land in the United Kingdom, usually receiving a greater proportion of inorganic P input versus manure containing P, has been demonstrated to have a greater ratio of P i :P o compared to improved grassland where a variety of P o forms exist (Stutter et al., 2015). Changes in agricultural land use from arable to grassland can potentially be used as a strategy to remediate soil and water pollution, with aims to reduce system P inputs and, in turn, losses, promote soil organic matter and sustainable P cycling amongst other benefits such as improving soil structure . Microbial cycling of P in soil by immobilisation or mineralisation is dependent upon the stoichiometry of C, N and P in the soil, those of the inputs to the soil, and how these relate to the stoichiometric requirements of the soil microbial biomass (Cleveland & Liptzin, 2007;Kuzyakov et al., 2019). Therefore, examining the molar ratios of C:N:P found in a soil over time may help to shed light upon the functioning state of the system of biogeochemical processes of nutrient cycling (Frossard et al., 2016).
The current study examines 50-to 70-year timescales of soil P change across multiple sites. Soil changes due to P fertiliser input changes have been studied at individual (or several) sites over decade to century timescales. Pasture receiving fertiliser for 13 years showed accumulated P o as broad phosphomonoesters groups related to biological P storage, but no specific change in myo-and scyllo-inositol hexakisphosphate, αand β-glycerophosphate and RNA mononucleotides (McLaren et al., 2017). A 30-year field trial in Switzerland comparing P fractions in more sustainable versus conventional farming techniques showed overall P i increases in the conventional farming treatment receiving mineral P fertiliser, but no change in P o using chemical fractionation or enzyme methods (Keller et al., 2012). When soils amended with differing forms of organic inputs for 62 years were compared using, 31 P NMR, Annaheim et al. (2015) found that no significant differences over time in P o quantities or forms but that P i had altered in amount. In Missouri, soil cultivation over 111 years showed an increase in resin-P and NaOH-extractable P i forms relative to native grass prairie on comparable soils, but opposingly, increased NaOH extractable P o under the grass (Motavalli & Miles, 2002). In summary, decades of fertiliser use commonly results in increased total P, often manifested in P i increases, but the added P forms may undergo conversions between P i and P o forms depending on soil conditions and potential for P o storage may become limited.
We compared archived topsoil samples from 35 agricultural sites in NE Scotland taken between 1951 and 1981 to those resampled from the same locations in 2017. Whilst this time comparison between samples has been generally compared in Lilly et al. (2020), we sought here to explore specific differences in functional forms of soil P. This allowed us to test the hypothesis that P forms (and by inference functions) have changed in agricultural soils over a period of 50-70 years, during which there have been substantial changes in UK agricultural practises. Both P i and P o change are explored using a range of extractants relevant to crop availability, leaching and organic P speciation. To better elucidate drivers and implications of change in P forms we explore soil properties related to P cycling, namely soil pH, organic matter content, macronutrient stoichiometry and P sorption strength. The available data preclude examining site-specific agronomic P balances or erosion and leaching components. Therefore, whilst we do not have detailed information about site-level spatiotemporal drivers of change, we are able to look across the sample population using broad grouping (start vs. end land use, soil properties) for P form changes, relationships between properties against the background context of agricultural change and the known mixed farming practises of the region, where soils are subjected to both arable and grazing over decadal timescales. This represents a long and important period of changing drivers acting on soil P change with potential to improve understanding future soil P trajectories.

| Sites and sampling
The samples used in this study (described in Lilly et al., 2020) were part of a resampling campaign of 37 experimental field sites located on working farms in NE Scotland, for which archived soil samples from earlier experiments were available. Their primary interest was change in soil carbon content over 50-70 years. A series of experiments were conducted between 1950 and early 1980s, by the Macaulay Institute of Soil Research (now The James Hutton Institute), to support national agricultural improvements in soil fertility and food security (Reith et al., 1987). Over 30 years, approximately 1000 field experiments were conducted on sites chosen to represent the NE Scotland soil conditions and responses to inorganic P inputs. Recently, 37 of these sites were relocated to the same sampling plot from original documents and resampled, producing a set of temporallypaired samples. The duration between the first original timepoint 1951-1981 and second resampling (T-RE; July-September 2017) was between 36 and 66 years. Climate change is not described at a site or regional level as this is beyond the present study; general change in the period is described in Figure 1.
The 35 locations used in the present study (summarised in Table 1) had arable cropping at the time of initial sampling. The most common regional system was a rotation of oats, root crop, oats/barley then grass. Root crops (swedes/potato) were used during the experimental phases as they were considered to give a detectable response to P inputs in the year of application. The average farm size was 20-32 ha and often incorporated an area of rough grazing. In contrast, the average farm size had increased to 71 ha by the time of resampling in 2017 and newer cropping habits favoured autumn sown arable crops such as oilseed rape and spring barley and a reduction in livestock. Table 1 groups land change between original and resampled times as (i) arable, remaining arable (n = 20), (ii) arable to current grassland (n = 12), where three sites were unclassified. Land change in the period intervening original and resampled times was not known and it would be expected that many fields experienced times of regular crop rotation or changed cropping with climate and market changes. Detailed site-specific attributes on P inputs, crop offtakes and P losses by runoff, leaching and erosion are not consistently available. Soil bulk density at resampling was available (Lilly et al., 2020) but not from original samples such that P indices were developed as mg P/kg basis not on mg P per soil volume. Agronomic categories stated in Table 1 use the soil test P results at original and resampling (but bulk density value at resampling for both timepoints) and show all sites at high or very high P status both times. Sites were dominantly freely drained, with four on imperfectly drained soil series, but the modification of field drainage is unknown. A range of NE Scotland's soil parent materials are reflected in Table 1  Soil P sorption class is a measure of the capacity of a soil to bind P as modelled from soil associations from 1, low to 2, moderate, to 3, high (SRuC, 2015). c The land use as noted from field survey according to: Sw, swedes; Tur, turnips; Pot, potatoes; SP Bar, spring barley; W Wheat, winter wheat; OSR, oil seed rape; Cer, cereals (used when stubble or ploughed land was indicated). This is then summarised in terms of two categories of change from arable to arable, or arable to grass, noting that all sites were originally arable and with no knowledge of change between the periods of survey.
d Agronomic soil P status compared for original and resampled timepoints as mg P/L volume of soil derived from modified Morgan's P concentration at timepoints and soil dry bulk density (using the value at time of resampling) according to SRuC (2015) as: high (H) = 13.5-30 mg P/L; very high (VH) = >30 mg P/L where categories of very low = <1.8 mg P/L; low = 1.8-4.4 mg P/L; moderate À = 4.5-9.4 mg P/L; moderate + = 9.5-13.4 mg P/L were not observed. The use of italics, plain font and bold indicates a decline, constant (no change >10%) and increase over time for a given site.
been used to develop a three-class system of P sorption capacity (PSC) used to modify the advisory soil P status per crop (SRuC, 2015); for example, a high PSC soil can attain a higher soil test P for some high P requiring crops. Sample distributions across the PSC classes were balanced (n = 10, 12, 13 for low, medium, high, respectively). Whilst we do not have site-specific soil properties (e.g. texture) the PSC system is derived from these and has a basis in contemporary agronomic P advice.
Archived T-OR samples from the original field trials were retrieved from the Scottish National Soil Archive and subsampled for analysis. These original dried samples were stored in dark, dry conditions, at room temperature. The resampling campaign for the T-RE samples, replicated the original sampling method and field locations of the T-OR samples. Composite samples totalling 1 kg per field site were collected using a screw augur in a W-pattern to a depth of 15 cm. These were stored in a cold room at 4 C before being air-dried at approximately 30 C and <2 mm sieved. Subsoils were not sampled but the general properties of the subsoils can be derived from the stated Association and PSC data in Table 1. All the topsoil samples (both T-OR and T-RE) were analysed simultaneously by current methods, to reduce errors associated with compiling old and new analytical method data.

| Soil extractions
Soil extractions targeting P i and P o across a gradient of readily soluble and labile, to refractory were conducted ( Table 2). Extraction of the soil samples with deionised water was conducted to determine soluble C, N and P forms, including orthophosphate which is considered the most readily plant available form of P (Sissingh, 1971). For each sample, 1 g of air-dried soil sieved to <2 mm was moistened with 1 ml of deionised water (Millipore) and left in the dark at room temperature for 22 h. A further 99 ml of deionised water was then added, and the samples were agitated in a reciprocal shaker (200 rpm, 1 h) before filtration (Whatman 42). Samples were stored at 4 C for <3 days before colorimetric analysis (Skalar San++) for dissolved organic carbon (DOC), nitrogen forms NH 4 + , NO 3 À , total dissolved N (TDN) and dissolved organic N (DON). Total dissolved P (TDP) and soluble reactive P (SRP) were determined by the same autoanalyser using the molybdate reaction, with and without an automated digestion procedure, allowing calculation of dissolved organic P (DOP) as the difference between TDP and SRP. Acid ammonium oxalate extraction was used to extract Al (Al ox ), Fe (Fe ox ) and P (P ox ) forms which are present in poorly crystalline forms or moderately bound to Fe-humus complexes, giving an indication of the soil's P sorption properties (Farmer et al., 1983). A Retch mixer mill with zirconia component was used to mill dry soil to a fine powder, which was then extracted (60 mg at 1 g:500 ml) in 0.2 M acid ammonium oxalate on a reciprocating shaker (200 rpm, 4 h in the dark at room temperature), then centrifuged at 5000g for 15 min. Extracts were filtered (0.45 μm) and stored at 4 C for <3 days before analysis by inductively coupled plasma optical emission spectroscopy (ICP-OES; Agilent 7500ce instrument; Agilent Technologies, CA, USA).
NaOH-EDTA extraction was used to determine cationchelated and organic matter complexed P in all samples with a subset carried forwards for solution 31 P NMR analysis (Cade-Menun & Liu, 2014;Cade-Menun & Preston, 1996). Samples of <2 mm, air-dried soil were extracted using 0.25 M NaOH with 0.05 M Na 2 EDTA at a soil to extractant ratio of 1:20 (wt/vol) at room temperature on a reciprocal shaker (200 rpm, 16 h) in the dark, centrifuged at 5000g for 15 min and then filtered (Whatman No. 42). An aliquot of each sample was diluted by 10 with deionised water and adjusted to pH 6.0-6.5 (1 M H 2 SO 4 ) for colorimetric analysis (as detailed above, denoted TDP NaOH-EDTA , SRP NaOH-EDTA , DOP NaOH-EDTA ) and analysis by ICP-OES (P NaOH-EDTA , Al NaOH-EDTA , Fe NaOH-EDTA ).
To quantify specific P, a representative subset of 10 samples (5 time-paired sites) was chosen for 31 P NMR analysis. A subsample of the NaOH-EDTA extract was lyophilised for each of these. To allow absolute quantification of P from the NMR analysis for this subset of 10 samples only, a milled fraction of the soil was analysed for total P using the NaOH-Fusion (P tot ) method (Smith & Bain, 1982).
In all three extractions, reagent blanks and an inhouse reference soil were included to quantify analytical error at appropriate stages in the processes. Previous analysis of the samples has provided data for soil pH (1:3 ratio of soil to water), exchangeable cations and MM-P (an extract of 0.5 M ammonium acetate/acetic acid, adjusted to pH 4.5 with colorimetric phosphate analysis) and for percentage elemental C and N concentration using flash combustion combined with mass spectrometry (Lilly et al., 2020). The MM-P method is the official agronomic P index used for Scotland's soils that are dominated by acid rock parent materials generally devoid of carbonates (MISR/SAC, 1985;SRuC, 2015).

| 31 P NMR spectroscopy of NaOH-EDTA extract solutions
Applying the advanced P speciation technique of 31 P NMR on a subset of samples (limited by resource constraints) allowed us to screen for major differences in specific P compounds with time that were not available via T A B L E 2 Description of phosphorus (P) extraction and chemical analysis methods used for all soils in relation to soil P functions, ordered in rows (top to bottom) from weak extractants (labile P, readily available pools) to strong extractants (including refractory pools)

Extraction and analysis method Analytical fraction defined Relation to P functions
Water extraction with colorimetric analysis including molybdate reactive P (Sissingh, 1971) Soluble reactive inorganic P (SRP) in mg P/kg Readily soluble P representing dominantly orthophosphate P i considered wholly and rapidly biologically available and/or available for leaching and runoff Total dissolved P (TDP) in mg/kg The total of dissolved P comprising SRP (dominantly P i ) and soluble organically complexed P o (by inference of being nonmolybdate reactive). This may be hydrolysed via enzymes for biological uptake and can indicate biological conversion of chemical fertilisers or organic P system inputs Modified Morgan extraction (MISR/SAC, 1985). An extract of 0.5 M ammonium acetate/acetic acid adjusted to pH 4.5. Operationally equivalent to the use of Olsen P method (0.5 M NaHCO 3 , pH 8.5) in less acidic soils Modified Morgan P (MM-P) in mg/kg. Agronomic P extraction for acid soils, government recommended method for Scotland A P fraction considered plant available to the next crop, including P both immediately available and the moderately bound P reserve in acid soil Acid ammonium oxalate extraction (Farmer et al., 1983). Strong acidic (pH 3) extraction for P, Al and Fe Oxalate P (P ox ) in mg P/kg Concentration of exchangeable P associated with amorphous surface binding complexes (principally Fe, Al) or bound to Fe-humus complexes P saturation (P sat ), calculated as a molar ratio where P sat = P ox /(Al ox + Fe ox ) The ratio of exchangeable P to the magnitude of surface exchange complexes that may be considered a P saturation proportional to intensity and amount of P capable of desorbing NaOH-EDTA extraction. Strong Alkaline extraction (pH 13.5). The addition of EDTA as a chelating agent greatly increases the extraction of organically complexed P and is commonly used before 31 P NMR analysis (Cade-Menun & Preston, 1996) Inorganic P by NaOH-EDTA (IP se ) mg P/kg Inorganic P from solubilisation of soil organic matter and removal of chelating cations complexing other forms indicating all weak to strongly bound inorganic P that is alkali soluble (likely excluding mineral bound P). May become biologically available directly or via biological exudates that counter complexing metal cations in soils Organic P by NaOH-EDTA (OP se ) as mg P/kg. Determined by difference as TP se À Ip se By inference of being non-molybdate reactive the organically complexed P that is alkali soluble. Longer-term P storagereplenishment pools. May become available in soils via biological complexing agents and/or enzymatic degradation pathways Total P by NaOH-EDTA (P se ) as mg P/kg Total P extracted by this method, as detected by inductively coupled plasma (ICP) quantification (TP se ) NaOH fusion total elemental digest and ICP quantification for total P of the solid phase sample (Smith & Bain, 1982) Total P as (P tot ) mg/kg Total P in the bulk soil comprising the sum of available, exchangeable and mineralbound P total and hydrolysable P measures on progressive extractant strengths. Samples were prepared for spectroscopy immediately before measurement by dissolving 150 mg of lyophilised sample into 1.5 ml of a solution of 1 M NaOH also containing 10% D 2 O for frequency locking and 0.2 mM methylene diphosphonic acid (MDP) as the internal standard for quantification of peaks. The solution containing the sample was centrifuged and then 1 ml was pipetted off avoiding solids which might contribute to line broadening and introduced into the NMR tube. Solution 31 P NMR spectra were acquired at 21 C on an Avance 500 II spectrometer (Bruker, Germany) operating at a frequency of 203 MHz using a 5-mm probe. Scans per sample numbered 512 with an acquisition time of 1.63 s. A recovery delay of 35 s was used. Spectra were processed in Bruker Topspin 4.0.8 software (Bruker Corporation, Billerica, MA). Line broadening of 10 Hz was used for all samples. Peaks were assigned as based on consideration of literature sources (Cade-Menun, 2005;Doolette et al., 2009;McLaren et al., 2015McLaren et al., , 2017Reusser et al., 2020;Schneider et al., 2016;Stutter et al., 2015), detailed in Table 3. Peaks were assigned at major group level: MDP standard (MDP, 16.8 ppm), phosphonates (P phon , 19-25 ppm), inorganic orthophosphate (P inorg , 5.3-6.3 ppm), phosphomonoesters (P mono , 3.2-5.3 ppm and 6.3-6.5 ppm), phosphodiesters (P di , 2 to À1 ppm) and pyrophosphate (P pyro , À4 to À5 ppm). We treated P mono group as a single integral, allowing standardised comparison within the study, but limiting detailed subdivision since we did not carry out spiking for precise isomer identification, deconvolution of confounded peaks Reusser et al., 2020), nor corrections for P di degradation products (such as α-glycerophosphate and β-glycerophosphate, a likely consequence of the alkaline extraction) occurring in the P mono region (Schneider et al., 2016). Quantification of peak areas (in mg P/kg soil dry matter) was determined relative to the quantified MDP internal standard and quantified as 89% (range 74%-108%). The total P determined by NaOH fusion method that was not accounted for by P NaOH-EDTA was designated alkali insoluble P (mean 26% in the NMR subset). Since the NaOH fusion total was carried out only on the 10 NMR samples, it must be noted that the TDP determined as P NaOH-EDTA on all samples represents approximately three-quarters of the total digest P.

| Data quality control and statistical analyses
The P saturation index (P sat ) was calculated as the molar ratio by P sat = P ox /(Al ox + Fe ox ). Mean pH values were calculated by log 10 [mean((1/10) ÀpH )]. Analytical data were checked for likelihood of systematic errors associated with indirect determination of DON and DOP that can result in considerable DON and DOP concentration uncertainty. We used a DIN:TDN versus DON diagnostic plot or SRP:TDP versus DOP diagnostic plot to assess this uncertainty. In this case, DIN = sum of NO 3 -N and NH 4 -N. This approach is based on Graeber et al. (2012), where it has been shown that at DIN:TDN <0.6 DON concentrations can be trusted; DIN:TDN ratios = 0.6-0.8 can show systematic deviations of the estimated DON concentrations from true DON; DIN:TDN ratios >0.8 indicate high random and systematic uncertainty; and DIN:TDN ratios >1 indicate overestimation of DIN, resulting in negative DON concentrations. We assume the same thresholds for DOP and SRP:TDN. Using this, we found all water extract DIN:TDN was <0.6 thus we have a high confidence in the use of DON data. In contrast, water extract SRP/TDP ratios of <0.6, 06-0.8 and >0.8 occurred in 31%, 29% and 40% of the original samples and 2%, 11% and 87% of those resampled. Because of this high likelihood of error (including 13 samples with SRP/TDP >1) we chose to reject further consideration of DOP in water extracts. The same checks found high confidence in using DOP data from the NaOH-EDTA extracts, whereby SRP/TDP of <0.6, 0.6 and >0.8 occurred in 77%, 14% and 9% of original samples and 43%, 54% and 3% of those resampled, respectively.
Statistical analysis was carried out in two ways: (i) general ANOVA analysis on all sites' data using Minitab (v16.2.4), (ii) multivariate analyses and analyses of the subset for 31 P NMR and were carried out in R (R Core Team, 2019;RStudio Team, 2019). The ANOVA (general linear model with Tukey differences) focused on the differences in timepoints (n = 2; original) with additional grouping factors for change in land use (n = 2) and soil PSC (n = 3; Table 1). These three grouping factors were approximately balanced (Table 1). No further records of liming, fertiliser, animal stocking records were available at the site-level with which to deduce locationspecific influences of agricultural change. Therefore, the changes between timepoints (and interactions of effects of land change and PSC on timepoint differences) were tested formally in the results, but then discussed as a whole population against the general context of background agricultural changes depicted in Figure 1. Pearson's correlation coefficient was calculated separately at both timepoints between all variables using function cor.test in R. The repeated measures correlation coefficient r rmcorr was calculated between all variables, across both timepoints. This method was used in addition to Pearson's correlation coefficient as it takes into account the nonindependence between the sites in repeated measures (Bakdash & Marusich, 2017). Principal components analysis (PCA) was conducted separately for each timepoint using function prcomp to investigate the main axes of variation in the dataset, and the variables associated with the loadings of each principal component (PC). To investigate the degree of separation between the timepoints linear discriminant analysis (LDA) was conducted. To reduce the effects of colinearity in the data, the LDA was conducted on the first six PCs of a PCA conducted on all data for both timepoints where site and timepoint were concatenated as a combined case, allowing subsequent calculation of the LDA with timepoint as the discriminatory factor. Data were log 10 transformed before testing according to Anderson-Darling tests and visual checks of quantile-quantile plots and histograms of the distribution of the data.

| General spatiotemporal context of the soils
The soil plots were relocated for resampling, but it remains important to understand the range in fundamental properties (spatial controls) and any differences between sampling periods (temporal control) as context for change in P pools and functions. The soil exchange complex defined by Al ox and Fe ox is an inherent property of the parent material and mineral weathering at longer timescales and was confirmed as not differing between timepoints (p > 0.05; Table 4). Very large ranges were observed that reflected the diversity of parent materials sampled by the study (for Al ox : 2416-13,241 mg/kg at T-OR and 1574-11,841 mg/kg at T-RE; for Fe ox : 3115-23,688 mg/kg for T-OR and 4067-14,940 mg/kg for T-RE; Table S1). The high significance of the PSC category term on Al ox , Fe ox in the three-way ANOVA (Table 4) relates to the fact that PSC is derived from average exchange properties modelled by parent materials, although no effects of interactions of PSC and time were observed. The total soil C (expecting inorganic C to be negligible in these acid parent materials) is an inherent property of soil landscape position and drainage with potential for management alteration and this was confirmed as not differing between timepoints (p < 0.05; Table 4). There was a threefold range in C content between sites consistently in time (2.2%-6.8% by mass for T-OR and 2.2%-7.5% by mass for T-RE) but no significance of the two broad land change groups in the three-way ANOVA (arable at T-OR, vs. either arable or grassland at T-RE) to suggest that land use influenced overall soil C. Soil total N is likely to be more influenced by management, such as the nature and amount of system N inputs and was slightly, but significantly greater in T-RE than T-OR (p < 0.05; Table 4). However, the soil total C:N molar ratio did not differ between timepoints (p > 0.05; mean ± 1 SE of 15.3 ± 1.9 at T-OR and 14.9 ± 6.4 at T-RE; data not shown). Hence, fundamental exchange complex and overall organic matter properties were consistent between timepoints at a spatial sample population level, not affected by the land change groups, but PSC was justified as a grouping factor. In contrast, soil pH, as a master soil indicator and control that is susceptible to management and other environmental change (e.g. climate), was highly significantly different between timepoints (p < 0.001; Table 4), being greater at T-RE (mean 5.9) than T-OR (mean 5.5), affected by PSC grouping (p < 0.01), with implications for control on labile P forms discussed below.

| Changes with time in analytically defined functional phosphorus pools
According to the explanation of P assays (Table 2) various individual techniques (not sequential extractions) progressively explored weakly to strongly bound P pools. The P pools readily soluble to water showed no significant differences in time for either the highly available SRP form or the sum of phosphate and complexed dissolved P defined by TDP (Table 4; p > 0.05). But TDP and SRP were significantly affected by the PSC grouping, where highest concentrations were associated with PSC 2 and lower in PSC 1 and 3. In contrast, the plantavailable MM-P decreased over time (p = 0.02, T-OR = 47 mg P/kg, T-RE = 37 mg P/kg).
Furthermore, no significant differences over time were observed in the P concentration of the soil exchange complex defined as P ox or in the P saturation of the exchange complex defined as P sat (Table 4, both p > 0.05). Whilst P ox and P sat were here derived from directly measured data the PSC was a modelled factor derived from mapped soil units. The higher PSC site groupings were supported by the greater directly determined values for P ox and P sat and vice versa. The NaOH-EDTA P extracted P by alkali-extracted organic matter and removal of chelating cations (Table 2; a mean 74% of total digest P on a subset, n = 10) found no difference over time in the TDP, SRP or DOP forms (Table 4; p > 0.05). The only difference in P forms with the land change factor was for P sat where the "Arable" group (arable at T-OR to arable at T-RE) of 0.090 was significantly smaller than the "Grassland" group (arable at T-OR to grassland at T-RE) of 0.105 (p < 0.05). The significant interaction between land change group and time showed that P sat varied significantly over time (T-OR and T-RE means of 0.075 and 0.106, respectively) within the land change group T A B L E 3 31 P NMR chemical shifts in ppm, for five sites for original and resample timepoints 1 (remaining arable) but did not differ in land change group 2 (means 0.104 and 0.105, respectively). The Pearson's correlation coefficient calculated at the two timepoints separately indicates strong positive correlation between SRP and MM-P in both cases (T-OR, r = 0.75 and p = <0.001; T-RE, r = 0.8 and p = <0.001), however, the repeated measures correlation coefficient (r rmcorr ) is much less positive and not significant (0.24, p = 0.16), reflecting more the complex nature of the relationship between these two variables in these data when including a temporal dimension, due to the opposite direction of changes in the means between the two timepoints (Table S2). P sat was calculated as positively correlated with SRP across both timepoints (r rmcorr = 0.81, p = <0.001) and less strongly but not significantly with MM-P (r rmcorr = 0.16, p = 0.36).

| Supporting data indicative of biogeochemical processes affecting P
A consistency in time in major soil factors of organic matter amount and exchange properties has already been presented, but extractions of macronutrients C and N and soil pH show strong differences over time. Using data from Lilly et al. (2020), soil pH (in water) was found to increase significantly between timepoints (p < 0.001, T-OR = 5.5, T-RE = 5.9). Significant greater values in T-OR than T-RE were observed in water extracted C and N measures, DOC (p < 0.001, T-OR = 1187 mgC/kg, T-RE = 212 mgC/kg), TDN (p < 0.001, T-OR = 75.9 mgN/ kg, T-RE = 40.3 mgN/kg), DON (p < 0.001, T-OR = 42.4 mgN/kg, T-RE = 27.6 mgN/kg) and NH 4 (p < 0.001, T-OR = 29.4 mgN/kg, T-RE = 6.7 mgN/kg). Water extracted NO 3 was at detection limits in many samples and not statistically tested. Hence, whilst water P form concentrations did not differ over time the macronutrient ratios governing biogeochemical processes, DOC:TDP and TDN:TDP, did differ significantly, showing that relatively P was less available than DOC and TDN in the soils at T-OR than T-RE. To investigate the C and N stoichiometry relating to organic P phosphorus, the molar ratio C tot :N tot :DOP NaOH-EDTA was calculated (T-OR = 330:21:1, T-RE = 354:25:1). The significant reduction in DOC between the timepoints was greater in magnitude than the relative reduction in concentrations of TDN. TDP increased relative to TDN. These changes are reflected in the change in molar ratio DOC:TDN:TDP (T-OR = 552:21:1, T-RE = 82:25:1) which indicates relatively higher concentrations of water-soluble nutrients in relation to DOC.
Ordination by PCA was conducted for the entire dataset of both timepoints (Figure 2), and for data for each timepoint separately. In the 'combined timepoint' PCA, 57% of the variance was described by PC1 and PC2 and the first four components described 79% of the variance. PC1 loadings are dominated at the negative end of the axis by a group of variables relating to soil pH and labile inorganic P forms and P saturation (pH, SRP, TDP, MM-P, P sat ) contrasting with a group at the positive end of the axis containing variables relating to Al, N and occluded forms of P (Al NaOH-EDTA , Al ox , TDN, N tot , P NaOH-EDTA, DON, P ox ). PC2 loadings show P-related variables grouping at the negative end of the axis (P sat , SRP, SRP NaOH-EDTA , TDP NaOH-EDTA ) contrasting with a group containing the water extracted measures of DOC, NH 4 and DON. LDA performed on the PC1-PC6 of this PCA, indicates separation between the timepoints. Using leave-oneout cross-validation, a misclassification rate of 3% was calculated for this model. Comparison of the PCAs conducted on the timepoints separately (data not shown) indicated that pH is an important factor explaining differences in the patterns of variance between the timepoints. The PCA for the original timepoint shows similar grouping of variables on PC1 and PC2 as described for the combined dataset. However, for the resample timepoint PC1 is dominated by pH alone, contrasting with a group containing N tot , TDN, and measures of occluded P (P NaOH-EDTA, P ox ). PC2 is dominated by a group of labile P forms (SRP, TDP, MM-P) contrasting with Fe ox .

| 31 P NMR spectroscopy for a subset of five sample time pairs
The 10 samples measured for NMR were similar in mean and range of P properties to the parent set, and thus representative. The subset mean and range values were: P NaOH-EDTA 990 mg/kg (578-1375 mg/kg); C tot 4.0% (3.1%-5.1%); P ox 1096 mg/kg (586-1609 mg/kg); Fe ox 8385 mg/kg (5209-18,198 mg/kg) and Al ox 5253 mg/kg (3119-9743 mg/kg). For SRP NaOH-EDTA , which was significantly greater in the resample timepoint data for the full dataset, all five sites showed greater values for T-RE than T-OR (T-OR mean 579 mg/kg, T-OR mean 794 mg/kg). Total P by the NaOH-Fusion (P tot ) method was measured in the range of 781-1833 mg/kg for the subset (Table 5; Figure 3). Four sites had greater values for T-RE than T-OR for total P by NaOH-Fusion (P tot ) and P NaOH-EDTA . The P NaOH-EDTA amounted to 44%-99% of P tot indicative of a varying pool of unextracted P (often termed alkaline insoluble P).
All 31 P NMR spectra ( Figure S1) were dominated by a peak at around 5.6 ppm designated as inorganic orthophosphate (Table 5). Only three samples showed detectable diester peaks (between 2 and À1 ppm) and seven samples pyrophosphate (À4 and À5 ppm). All samples showed significant combined patterns of peaks between 5.3 and 3.2 ppm, identified as the phosphomonoester region. For our screening level of interpretation peaks in this region were identified (depicted Figure S1, but not confirmed by spiking so not quantified) as various inositol-hexaphosphate (IHP) isomers including myo-IHP and scyllo-IHP. Peaks for diester degradation products were also tentatively identified as α-glycero-phosphate ($4.6 ppm) and β-glycerophosphate ($4.2 ppm) although these were often not obvious sharp features in these spectra. No samples showed obvious phosphonate peaks (expected between 19 and 25 ppm if present). All sites showed an increase over time in P inorg quantified by 31 P NMR, with increases ranging between 55 mg P/kg to 714 mg P/kg. In contrast, the differences over time in the P mono signal were mixed with two sites showing decreases (À58 and À563 mg/kg) and three sites showing increases (35, 59 and 150 mg/kg). A pyrophosphate peak was identified in two of the original timepoint samples (8-10 mg/kg) but all five resampled (between 2 and 9 mg/kg). In summary, screening a subset of samples by 31 P NMR showed that similar P compound diversity occurred in both original and resample timepoints, except for pyrophosphate which occurred at detectable levels in all resample spectra but was not observed in original samples for three of the sites and the consistent increase in P inorg .

| DISCUSSION
We investigated whether quantities and forms of P had changed in agricultural topsoils that were sampled at two timepoints 50-70 years apart having been relocated at field and plot level. This period represented one of substantial change in regional farming practises (Table 1) and potential mechanisms of change affecting the P forms are discussed in relation to this agricultural change context and supporting data on soil macronutrients and pH.
First, we were able to establish that the two timepoint populations of soil samples were broadly similar in their dominant soil properties associated with organic matter amount (C content) and the Fe and Al exchange properties denoted by Fe ox and Al ox . This gave confidence that we were comparing over time the same soil population in terms of major soil formation factors associated with landscape position and parent materials (operating at hundreds of years timescales). The soils were all originally, and remain at, high soil P level (Table 1) and can be considered to be P enriched with large historical P inputs. The main findings were that few of the P extract concentration measures here showed changes in time across the whole population (Table 4). No change in time was found for water extract P (highly mobile), for oxalate extracted P (P in the geochemical exchange complex) and for the P forms determined using the NaOH-EDTA extraction (P released from organic matter). The MM-P measure, used as the national soil test P index in Scotland (SRuC, 2015) did decrease over the 50-70 year gap, which contrasts with a more recent analysis by Edwards et al. (2016) who found no trend in MM-P between 1993 and 2010. There was a large and significant decrease in water-soluble DOC and N concentrations over time that contributed to higher TDP relative to DOC and TDN in contemporary samples. A significant increase in soil pH dominated the biplot differences between T-OR and T-RE (Figure 2), illustrating the influence of this master variable of biogeochemical processes on overall patterns of change across the dataset, and by inference on the functioning of the system of nutrient cycling in the soil.
Although all sites have remained in mixed agricultural use throughout the intervening time between sampling, styles of farming and farm structure have changed between the measurement points used in this study (Figure 1). Since the 1990s fertiliser recommendations in Scotland, as elsewhere, have sought a balance between fertiliser inputs and crop off-take, resulting in corresponding drops in rates of P application, which is reflected in an overall trend in reduction of the agronomic measure (MM-P), between the T-OR and T-RE samples (Lilly et al., 2020). The change in fertiliser practises manifested themselves in the study region as decreases in P fertiliser application rates of up to half for turnips, spring cereals and grassland, but steady inputs to potatoes and increased usage on the newly appearing winter cereals and a steady contribution of manure to the overall P inputs (Domburg et al., 1998). This has been a response to continued regulation at European to study region level of the quantity and timing of agricultural P inputs, both chemical fertilisers and efficient usage of manures, to reduce system P excesses and avoid ecologically damaging P loss (Crooks et al., 2019). Aside from the decrease over time in MM-P, it is surprising that the water soluble, P ox and NaOH-EDTA extracted P concentrations have shown no differences. A statistically significant positive correlation between SRP and MM-P within timepoint was found in the current study. The observed decline in MM-P (decreasing T-OR to T-RE) contrasting with no significant change over time in water soluble SRP may relate to pH effects of the extractions. The increase in pH of the soils with time may have increased SRP solubility in unbuffered water extracts masking any decline in water soluble P at a fixed pH, as was detected with the MM-P extraction that is buffered at pH 4.5 . Despite overall regional decline T A B L E 4 Phosphorus fractions and supporting data mean and variation (1 SE) for original and resampled soils including three-way ANOVA testing results with group factors of time (n = 2), land change (n = 2) and soil association derived P sorption capacity (n  in areas of grassland favouring land use change to cereals (Figure 1), the sample population included sites that were originally arable and were grassland on resampling. However, this simplistic (two-point, without knowledge of intervening period) analysis of land change as a factor showed negligible impact on P pools. The exception of the smaller P sat at the original timepoint in sites that remained arable (then a subsequent increase in time), and larger P sat at original timepoint in sites that latterly changed to grassland (no resulting change in time), cannot readily be explained without further site-specific detail.
Limited change in P may be explained by the slow hysteresis in the trajectory of declines in different functional P pools relative to P input declines that is commonly associated with the 'legacy P effect' (Rowe et al., 2016). Whilst mineral bound P is subject to longer-term mineral weathering release, other moderately refractory components of total P (e.g. as defined by NaOH-EDTA total P) act to replenish labile P pools (e.g. 'intensity' and 'resupply' in water soluble P or MM-P) . Thus, reversal of historical P loading requires depleting labile pools to draw-down total P. Limited differences in time across P extractant strengths demonstrate the long timescales involved in recovery for high P status agricultural soils. Within arable rotations periods of potatoes and winter cereals, both high fertiliser P inputs, will also have maintained high contemporary P inputs. In the study region, the lag effect between changes in P fertiliser recommendations and response of soil P status has previously been noted . Soil pH latterly is also likely to have strong influence on P mobility as pH increases P solubility (the 'intensity' of the labile pools; Lumsdon et al., 2016). Liming to bring acid soils to pH 5.5-6.0 favouring growth of crops, especially cereals and balance trace element availability, was a subsidised activity until the 1980s. We do not have evidence for change in regional liming practise, nor specific data for these plots. The increased current soil P contradicts the earlier liming subsidy removal, but it is also possible that soil testing, advice and general awareness is also greater currently.
Quantification on specific P compounds on a subset of samples by 31 P NMR showed increases in the amount of P inorg for four out of five sites, and five of the T-RE samples showing inorganic pyrophosphate peaks while these were only noted in two of the T-OR samples. Results for the organic P mono region were more mixed with three sites with greater results for this region in T-RE than T-OR, and two with lesser values. Soils subject to agricultural amendments of a range of types (chemical fertilisers and manures) often show the greatest P accumulations comprising inorganic ortho-phosphate, but this may be cycled to organic P storage forms often considered to end up in the more recalcitrant monoester P pools (Annaheim et al., 2015;Keller et al., 2012). However, there is evidence that a maximum upper limit exists for P stored in organic forms when P is administered in excess of the agronomic optimum (McLaren et al., 2017). The mixed results for the changes in the P mono region could potentially indicate that the maximum limit for storage of P in organic forms has been reached for these agricultural soils which have been in consistent longterm use. The screening of population changes in time using a limited set of NMR analyses here showed potential for the technique to look beyond the operationally defined pools associated with total P by selective extractants and reactant combinations. Further work should target mechanisms of change led by site-specific data on changes in inputs such as manure versus fertiliser, lime and field management such as cultivation.
Considering macronutrient stoichiometry, we found C tot :N tot :P tot molar ratios (T-OR = 85:6:1, T-RE = 76:5:1), for the subset of 10 samples used for 31 P NMR analysis (where total digest P was available) that were similar to previously found for agricultural soils. A measure of total digest P was not available for the dataset as a whole, however we calculated molar ratios of C tot :N tot :DOP NaOH-EDTA of our samples (T-OR = 330:21:1 and T-RE = 354:25:1) and compared these to published data. In a set of three studies over 60 years in West Africa, Frossard et al. (2016) reported differences relative to a control (C tot :N tot :DOP NaOH-EDTA 320:33:1) of 181:18:1 (low application of mineral fertiliser), 396:38:1 (low application of mineral fertiliser and low manure application), 234:22:1 (high application of mineral fertiliser) and 333:32:1 (high mineral fertiliser and high manure). This comparison indicates that the nutrient molar ratios found in our study resemble other long term cropped systems with high inputs of manure and mineral fertiliser, and that differences in C tot :N tot :DOP NaOH-EDTA molar ratios between timepoints may be functionally negligible, in terms of overall amounts of C, N and P present. However, we found differences in labile water extracted measures indicating possible differences in nutrient cycling conditions at the two timepoints. Substantial decreases in water extractable DOC and DON in the T-RE samples than T-OR (as opposed to total elemental compositions) indicate change in available forms of organic matter, and likely how it is cycled by the soil microbial population. The decrease in water soluble N over time corresponds with declines in fertiliser N inputs generally in the region (Figure 1) but contradict a significant but small increase over time in soil total N (Table 4). For soil total C no change over time observed in the present study corroborated with other findings (Buckingham et al., 2014). The lowering of DOC:TDP and of TDN:TDP ratios over time (Table 4) indicate that relative to C and N, the P available to microbes was more depleted in T-OR than T-RE. According to Stutter et al. (2018), the DOC:TDP ratios at T-OR are elevated more like semi-natural soils and within a range of balanced resource to heterotrophic microbes (C:P $50 to $900), whereas those at T-RE depict more current global agricultural soils enriched in N and P. Amongst the effects related to this could be that microbes cycled the soluble P in the soils during 50-70 years ago more tightly in the presence of excess C and N, whereas in contemporary soils the P is more in excess to C and N. As such P in excess may be stored in soils in more P i than P o forms (Stutter et al., 2015), this may partly explain the accumulation of orthophosphate P i when measured specifically with 31 P NMR, although the observation was limited to a subset of soils. A further implication of changing DOC:TDP ratios is that dissociated DOC competes with phosphate sorption on anionic exchange sites. Lower ratios at T-RE could favour stronger retention of P on soil surfaces and have contributed to the greater P sat for resampled soils within the subset of those which were arable at both T-OR and T-RE. Note: For chemical shifts used see Table 3. Abbreviations: ICP-OES, inductively coupled plasma optical emission spectroscopy; nd, below detection limit; P, phosphorus.
A possible confounding issue in these P results is that an unknown degree of difference could be attributed to alteration of the T-OR samples during dried state in the soil archive, despite that this practise is overall considered an effective protocol to give soil resources for time comparisons (Blake et al., 2000). Air-dried soils stored for 3 years were observed to have decreases in inorganic P and increases in organic P as extracted by NaHCO 3 , these changes being attributed to disruption of organic matter on mineral surfaces, solid phase diffusion of phosphate into soil particles and microbial cell decomposition (Turner, 2005). These results suggest that the T-OR soils might be expected to be show greater values for labile forms of inorganic P than the T-RE samples. However, no change in P forms but considerable differences in DOC and dissolved N forms were observed. It is important to note that a storage effect on the pH of the samples cannot be ruled out. Previous studies have shown that the direction and magnitude of any such effects, if present, are not easily predictable (Prodromou & Pavlatou-Ve, 1998). Another factor is the effect of rewetting which has been shown to cause changes in extracted P, some of which is due to microbial activity (Blackwell et al., 2009). Although the soils at both timepoints were dried using a similar method and were extracted together for all measures identically, those which had been stored for less time may have exhibited different microbial reactions to rewetting with an unconstrained possible effect on P forms extracted. This may also have affected the greater concentrations of DOC and water extracted N in the T-OR relative to T-RE samples. It may have been possible that some microbial activity continued in microsites during earlier years of the dried soil storage. Drying and rewetting is known to affect microbial behaviour relating to both C and N and it is possible that this too had a different effect in the T-OR versus the T-RE samples, for example, if the microbial community was more active in the T-RE samples it may have sequestered the DOC as an easily available energy source during the incubation period of the water extraction method, before measurement. Long term dried soil samples represent an important resource of information on soil states over time, so further studies investigating the nature of storage effects, especially on more labile nutrient fractions, are recommended so allow full benefits of this resource.

| CONCLUSIONS
The comparison of topsoils between 1951 and 1981 to 2017, allowed testing of the hypothesis that P forms (and by inference functions) have changed in high soil P status agricultural soils over 50-70 years of substantial changes in agricultural practises. Our comparison of soils over time used basic grouping factors (land change and soil PSC) but without individual site fertiliser and cultivation management histories. A significant decline over time in the national agronomic P test result was linked to documented declines in chemical P fertiliser for the United Kingdom and the study region and steady manure P contributions. Consistent with legacy P concepts, a slow depletion of labile and moderately refractory P pools (that replenish labile P) seems to be occurring in response to the system fertiliser input declines that only happened latterly in the studied period. This occurred against a background of increased soil pH between the original and resampling F I G U R E 3 Mean 31 P NMR quantification results in mg/kg dry matter, for five sites at original (black bar) and resample (grey bar) timepoints. Error bars are SD of the mean timepoints of 0.4 pH units, likely also due to altered management. The result of no significant overall change in time observed for water soluble P, exchangeable P (by acid ammonium oxalate), nor complexed P (by NaOH-EDTA) was surprising given the multidecadal duration and that at least 12 of the sites experienced change from arable to grassland. Screening a subset of five sample pairs for specific P compounds by 31 P NMR showed no consistent change in time for diversity of P compounds except for a consistent increase in orthophosphate P in contemporary samples. In addition, land use changes associated with mixed land use practises appeared not to have significant effects on P, although our data are limited by specific field data especially on intervening period. This has implications for future land management seeking remediation of high P legacy and prompts further work across timescales in mixed agricultural use systems, including monitoring of P using additional measures to the traditional agronomically relevant P index. This study represents a long and important period of changing drivers acting on soil P change with potential to improve understanding future soil P trajectories. The value of archived historical soils samples is shown, if care in interpretation of storage implications is made.