Phosphorus buffering determines how soil properties and rainfall influence wheat (Triticum aestivum) yield response to phosphorus fertiliser

Field experiments

Forty P rate-response field experiments were completed with the same experimental design. The design included five rates of P fertiliser: 0, 5, 10, 20, and 40 kg P ha⁻¹ (referred to as P0, P5, P10, P20, and P40, respectively) drilled at about 3 cm below the seed as mono-ammonium phosphate (MAP). The field experiments were arranged in randomised complete block design with three or four replicates. The rate of nitrogen (N) applied as urea at sowing was adjusted for each P rate to ensure all plots received 20 kg N ha⁻¹ as fertiliser at sowing. The appropriate amounts of urea (46% N) and MAP were drilled below the seed at sowing together with 35 kg ha⁻¹ potassium sulfate (15 kg potassium (K) ha⁻¹, 6 kg sulfur (S) ha⁻¹).

Agronomic management of the field experiments was consistent with current cropping systems in WA. Wheat (Triticum aestivum L.) cv. Scepter, a commonly grown variety in WA (Shackley et al. 2021), was sown at all sites. Experiments were sown with plot seeders fitted with no-till tines and press wheels. The experiments were sown in fields used for long-term no-till crop production, where the preceding crops were wheat (16 experiments), canola (Brassica napus L.) (13 experiments), barley (Hordeum vulgare L.) (three experiments), lupin (Lupinus angustifolius) (three experiments), pasture (three experiments), and oats (Avena sativa) (two experiments). Only one crop per year is grown in this region. Weed control was achieved by applying mixtures of knockdown and pre-emergent herbicides immediately before sowing (e.g. Harries et al. 2020) and selective herbicides post-sowing as required. The median and interquartile range for sowing date was 18 May, and 11−27 May, respectively, which fall within the window used for sowing wheat by grain producers in WA (Fletcher et al. 2016; Harries et al. 2020). Post-sowing N was surface-applied as urea or urea-ammonium-nitrate at a rate appropriate for each site based on assessments of yield potential and soil N supply. The median and interquartile range for total N application was 72 kg N ha⁻¹ and 62 to 96 kg N ha⁻¹, respectively, similar to the rates regularly applied to cereals in WA (Planfarm 2021). The size of plots sown differed slightly between the two providers managing the experiments. Half of the experiments were sown by CSBP Field Research (https://csbpfertilisers.com/advice/research-and-agronomy/field-research) and plots were 15 m long and 1.89 m wide, at 0.27 m row spacing, and half were sown by Summit Field Research (https://www.summitfertz.com.au/research-agronomy) and were 10 m long and 2.25 m wide at 0.22-m row spacing.

Soil measurements

Soil at the sites were classified according to the Australian Soil Classification (ASC) (Isbell R.F. and the National Committee on Soil and Terrain 2021). The ASC soil orders at the field experiment sites were Tenosols (15 sites), Chromosols (11 sites), Kandosols (11 sites), Calcarosols (two sites), and Kurosols (one site). For a brief description of the soil orders, see Supplementary Table S1.

Soil samples were collected from 0–5, 5–10, 10–15, 15–20, 20–30, 30–40, and 40–60 cm depth by combining five cores at each depth from plots with no fertiliser added in each replicate. Soil samples were dried at 50°C in a forced draught oven and sieved (<2 mm) for soil physical or chemical measurements.

Soil chemical properties were analysed according to Rayment and Lyons (2011). Soil pH (CaCl₂) (Method 4A1), Colwell extractable P (Method 9B), P buffering index (PBI) (Method 912C), Colwell extractable K (Method 18A1), KCl-40-extractable S (Method 10D1), oxalate-extractable iron (Ox-Fe) and aluminium (Ox-Al) (Method 13A1), electrical conductivity (EC) (Method 3A1), and organic carbon (OC) (Method 6A1) were measured. Plant-available P using the diffusive gradient in thin films method (DGT-P) (Mason et al. 2010) was also measured.

Particle size analysis of the fine earth fraction (<2 mm) was measured using the pipette method (Indorante et al. 1990). Soil texture was determined using the soil texture wizard package in R (Moeys 2018). Particle size distribution for gravel (≥2 mm) was measured using FAO size classifications (Jahn et al. 2006).

Rainfall measurement

Rainfall data were recorded for each site. At each site, a logging pluviometer (ICT International, Armidale, NSW, Australia) was installed immediately after sowing, and daily rainfall was calculated from each recorded 0.2 mm tip and its associated time period. Complete annual rainfall data for each site were created by merging observed rainfall and the interpolated data from the nearest rainfall station in the SILO database (Jeffrey et al. 2001).

Data transformation and statistical analysis

Data transformation, analysis and visualisation was done in R (R Core Team 2021). Imputation of left-censored values was done using Iterative Robust Model-based Imputation (Kowarik and Templ 2016). Left censored values are those below the reporting threshold for the laboratory used. The threshold used, and percent of samples in the soil data below this were: Colwell P (2 mg kg⁻¹, 4%); Colwell K (15 mg kg⁻¹, 2%); and DGT-P (4 μg L⁻¹, 10%). All the left-censored data were from soil samples below 10 cm depth. Linear and quadratic models were fitted to grouped data with lmList function (Bates et al. 2015). The arcsine-log calibration curve (ALCC) was fitted using the soiltestcorr package (Correndo et al. 2023). Figures were produced using the ggplot2 and ggthemes packages (Wickham 2016; Arnold et al. 2021) in R.

Relative yield (RY) was calculated for each site using the mean values of grain yield for the P treatments. Analysis of variance and comparison of treatment means using least significant difference (5%) showed grain yield was lower at P40 than P20 at five of the 40 sites. To avoid a reduction in yield at P40 compared to P20 leading to an over-estimate of relative yield, the maximum yield from the P treatments was used to calculate RY (%):

Mean values for soil chemical or physical properties were calculated to assess the relationships between crop response and soil test data from different depth intervals (e.g. simulated sampling depths of 0–5, 0–10, 0–20 cm). The depth-weighted mean was calculated based on the depth of the sampling intervals (e.g. Sharpley 2003). The values for soil chemical properties were not adjusted for soil gravel content.

Description of the database

The geographical distribution of the field experiments is in Fig. 1. The sites provided a good geographical representation of the agroecological zones in WA. Median growing season rainfall for the experiment sites across 4 years was 223 mm, and the range was 135 to 496 mm (Fig. S1). The sites also captured the north-south gradient in mean annual temperature (ranging from 15 to 21°C) that is typical for this region (Harries et al. 2021). Summary statistics for site mean grain yield, RY and the absolute yield response (max grain yield – grain yield at P0) are in Fig. S1.

Fig. 1.

Geographical distribution of the field experiments in the grain-producing region of Western Australia (shaded grey) from 2018 to 2021. Rainfall isohyets are in mm per year, based on rainfall observations from 1988 to 2017 (BoM, Bureau of Meteorology). Trial sites are labelled based on the closest town.

Summary statistics for selected soil properties are in Fig. S2. Median clay content increased as soil depth increased; from 9 to 22% at 0–5 cm and 40–60 cm, respectively. Similarly, median gravel content increased with soil depth, from 4% at 0–5 cm to 20% of soil weight at 40–60 cm. Most of the soil profiles were acidic, with a depth-profile typical for subsurface acidity induced by crop production; pH_Ca at 10–20 cm was lower than the soil above and below this depth (e.g. Tang and Rengel 2003). Median soil pH_Ca at 0–5, 10–15, and 40–60 cm was 5.6, 4.9, and 5.7, respectively.

The distribution of PBI values for 0–10 cm in this study was positively skewed, with 38% of values in the Moody (2007) ‘very very low’ category of 15–35. The experimental sites used in this study were representative of the PBI values that occur in this region; the frequency distribution of PBI values is like that from a large-scale survey (Weaver and Wong 2011) of 109,000 samples from 0 to 10 cm in the grain-producing region of WA (Fig. S3). The PBI showed an increasing trend with depth, whereas Ox-Al was more uniformly distributed with depth. The OC, Colwell P and DGT-P were stratified, with the highest values in the surface 10 cm (Fig. S2).

Our preliminary analysis showed that relative yield was not related to sowing date, total nitrogen fertiliser applied (kg N ha⁻¹) and preceding crop, and these variables were not included in the data analysis reported here. Relationships between these variables are in Fig. S4.

Approach to statistical analysis

The steps taken to address our research objectives are described in detail below. In some cases, the results from one step in our analysis was used to guide the analysis in the following step.

Objective 1. Assessment of soil properties that had the greatest influence on relative yield

Objective 2. Investigation of how climatic conditions influence yield response to P fertiliser

The effect of rainfall on the yield response to P fertiliser was assessed using the sliding window approach (Bailey and van de Pol 2016). The purpose of this approach is to examine whether there are climatic conditions within time-based windows that relate to an observed response. In this study, we analysed relative yield in relation to the total rainfall within time-based windows using a function written in R. The opening time for the windows ranged from −20 weeks to 19 weeks relative to sowing, and the closing time ranged from the opening time from +1 week to 20 weeks after sowing. For example, for the opening time of −2 weeks relative to sowing, windows closing from −1 week to 20 weeks relative to sowing were assessed. The sliding window analysis included fitting linear and quadratic models to all RY-rainfall data, or groups of data based on the PBI threshold identified above (PBI 0–10 cm <56 (PBI <56), PBI 0–10 cm >56 (PBI >56)) for each time-based window. Akaike information criterion (AIC) was used to compare different model fits because it accounts for goodness of fit and the number of model parameters. The model with the lowest AIC is the preferred model.

Objective 3. Examine whether there were interactive effects between rainfall and soil properties on relative yield

Interactive effects between rainfall and soil properties on relative yield were assessed using regression tree modelling. The modelling was done after the analysis described in Objective 2 was completed, where a systematic relationship between rainfall before sowing and relative yield was identified for sites with PBI was <56.

Objective 1. Assessment of the soil and site properties that had the greatest influence on relative yield

Investigating whether soil chemical profiles differ between sites clustered by relative yield

ANOVA of soil properties at each depth revealed that differences between RY groups mostly occurred in the surface 20 cm (Fig. 2). The lowest RY group (Group 1) was associated with the highest PBI, Ox-Al and OC values, the lowest DGT-P and the ratio Colwell P:PBI. The only instance where soil properties for RY Groups 2 and 3 differed was for percent gravel in 19–63 mm fraction at 5 to 10 cm; it was lower in RY group 3 (0%) than Groups 1 and 2 (0.85 and 0.54%, respectively).

Fig. 2.

Soil chemical profiles for relative yield clusters. Depths where soil properties differed significantly (P < 0.05) between relative yield groups within a depth layer are marked with *. PBI, phosphorus buffering index; Ox-Al, oxalate-extractable Al (mg kg⁻¹); OC, organic carbon (g kg⁻¹). Units for Colwell P and DGT-P are mg kg⁻¹ and μg L⁻¹, respectively. Pct. 19–63, percent gravel (19–63 mm) of whole soil by weight.

Investigating soil test calibration curves

An examination of the relationships between RY and soil chemical or physical properties showed closer relationships between soil properties from combined depths than the individual depths. A quadratic model was fitted to all combinations of soil depth and soil property data (not shown here). The better model fits were observed for soil data 0–10 cm for DGT-P, OC or Ox-Extr. Al and relative yield with R² ranging from 0.19 to 0.22. However, RY was not closely related to any single soil property, so we focussed on combinations of PBI and measurements of plant-available P (Colwell P, DGT-P).

To investigate whether PBI influenced the relationships between RY and Colwell P or DGT-P, the data were analysed based on PBI classes developed for Australian soils by Moody (2007) using the ALCC approach. Our initial analysis (data not shown) indicated insufficient data points (<8) in some of these PBI classes to fit a calibration curve, so the six PBI classes proposed by Moody (2007) were aggregated into two classes; i.e. 0–35 (n = 22) and >35 (n = 18). The median r-values (and 90% CI) for the calibration curve for Colwell P were low for the whole data set (r=−0.07, −0.35–0.2), PBI 0–35 (r = 0.02, −0.39–0.41), and PBI >35 (r=−0.11, −0.43–0.26). The r-values for the calibration curve for DGT-P was highest for PBI >35 (r = 0.63, 0.52–0.83), followed by the whole data set (r = 0.45, 0.23–0.63) and PBI 0–35 (r = 0.34, −0.10–0.68). The PBI classes presented by Moody (2007) were for a national data set and may not be optimal for WA soils; therefore, we investigated how the r-values changed for different PBI ranges.

Overall, the goodness of fit for soil test calibration curves was best for DGT-P at the high end of the range of PBI values observed in this study (Figs S5 and S6). For DGT-P, the best calibration was for a sampling depth of 0–10 cm and PBI range 56–396, and the median r-value was 0.95 (90% CI was 0.86–0.98) (Fig. 3). There were seven other combinations of sampling depth and PBI range with a median r-value >0.9, but all had PBI 0–10 cm greater than 56 and sampling depth within 30 cm depth. The goodness of fit of the calibration curve for DGT-P decreased as the range of PBI values increased; for 0–10 cm sampling depth and PBI range 45–396, median r-value was 0.83 (0.65–0.97), and PBI range 40–396 the median r-value was 0.72 (0.48–0.95) and this trend continued as the lower limit of the range of PBI decreased further (Fig. S6). The 90% CI for r-value for a lower limit for PBI of 37 or below did not overlap with that for PBI 56–396. In this dataset, there appears to be a threshold of PBI 56 for DGT-P; relative yield at sites with a PBI below this threshold were not as closely related to DGT-P as those above. The best calibration curve for Colwell P was for a sampling depth of 0–30 cm and a PBI range 9–17 (n = 10), and the median r-value was 0.64 (0.04–0.9) (Fig. 3). There were several calibration curves for Colwell P that had similar r-values as the best calibration curve, but all were for a narrow range of PBI (e.g. 9–17, 9–18 or 9–19) but included a wide range of sampling depths (from 0–20 cm to 0–60 cm). In comparison to DGT-P, the 90% CI for r-values was wider for Colwell P that for DGT-P. Based on the overlap of the 90% CI for r-values for Colwell P for different sampling depths and PBI ranges, we could not clearly identify a best sampling strategy for Colwell P. It is worth noting that the best calibration curves for DGT-P and Colwell P only accounted for 20 sites, and that RY was not closely related to Colwell P or DGT-P at the other 20 sites.

Fig. 3.

Soil test calibration curves with highest r values identified for Colwell P (a) and DGT-P (b) from the analysis of soil test calibration curves for differing PBI ranges and soil sampling depth. The highest r-value for DGT-P was found for a sampling depth of 0–10 cm and PBI 0–10 cm range 56 to 396. The highest r-value for Colwell P was for a sampling depth of 0–30 cm and PBI 0–10 cm range 9 to 17. The r-values and critical soil test value (CSTV) (for 90% relative yield) were calculated using the arcsine-log calibration curve method. The median and 90% confidence interval (CI) for the calculated parameters are for the bootstrap samples done for these calibration curves.

Objective 2. Investigation of how climatic conditions influence yield response to P fertiliser

The sliding window approach did identify relationships between rainfall and RY, but only when the sites were analysed using the PBI groups identified above. The quadratic models provided a better fit than the linear models, and was used to identify time-based windows with the closest relationship between rainfall and RY. For sites where PBI <56, the best fit was for the window from −5 weeks to −2 weeks relative to sowing (Fig. 4a). For PBI >56, the best fit was for the window from −8 weeks to 12 weeks relative to sowing (Fig. 4b). A matrix plot of AIC values for different combinations of opening and closing times (Fig. S7) shows the best model fits for PBI <56 were for time-based windows that opened in the weeks just before or after sowing, providing some evidence that timing of rainfall is important for RY. Using a threshold of delta AIC ≤5, there were six windows with a similar model fit to the RY-rainfall data from −5 weeks to −2 weeks; the opening values ranged from −5 weeks to −7 weeks, and the closing values ranged from −1 week to 8 weeks relative to sowing (Fig. S7a). For PBI >56, the best fit was for the window from −8 weeks to 12 weeks relative to sowing (Fig. 4b), and there were 25 windows with a similar model fit as weeks from −8 weeks to 12 weeks; the opening values ranged from −13 week to 0 week, and closing values ranged from 3 weeks to 13 weeks relative to sowing (Fig. S7b). For PBI >56, the effect of timing of rainfall was not apparent.

Fig. 4.

Relative yield–rainfall relationship for time-based windows: (a) from −5 weeks to −2 weeks relative to sowing for sites where phosphorus buffering index (PBI) 0–10 cm <56; or (b) from −8 weeks to 12 weeks relative to sowing for sites where PBI 0–10 cm >56. AIC, Akaike information criterion.

Objective 3. Examine whether there were interactive effects between rainfall, soil properties and relative yield

The data for the window from −5 weeks to −2 weeks relative to sowing were analysed further to examine whether soil properties influenced the RY–rainfall relationship. Because our regression tree analysis showed soil properties at 0–10 cm were the most important for predicting RY, we used the data for this depth only. Regression tree analysis was done to examine whether there were interactions between soil properties at 0–10 cm and rainfall from −5 weeks to −2 weeks relative to sowing. The first split was rainfall and the second was pH_Ca, and the resulting R² was 0.59 (Fig. 5). There was no relationship between rain in this window and Colwell P 0–10 cm, or Colwell P 0–30 cm (data not shown).

Fig. 5.

Regression tree for relative yield, based on rain in a time-based window from −5 weeks to −2 weeks relative to sowing for sites where phosphorus buffering index 0–10 cm was <56, and on soil properties 0–10 cm. Values in coloured boxes are the mean values for relative yield and the percentage of observations that fall within that node.

Our analysis of rain, soil and relative yield data for the window from −8 weeks to 12 weeks relative to sowing for the PBI >56 data showed a relationship between rain in this window and DGT-P. The DGT-P decreased from approximately from 50 to 10 μg L⁻¹ as rain during this period increased from 150 mm to 350 mm (data not shown). However, it is unlikely that this is a causal relationship because the window closed from 2 to 3 months after the soil sampling was done, meaning that some, or all the rain occurred after DGT-P was measured.

Objective 1. Assessment of the soil properties that had the greatest influence on relative yield

The PBI threshold we identified for determining whether rainfall and soil pH_Ca, or DGT-P should be used to predict yield response to P ought to be viewed as an indicator, rather than a precise cut-off point. The PBI threshold of 56 was determined empirically by a search for the soil test calibration curve with the highest r-values over a range of soil sampling depths and PBI values in the current study. However, previous work in WA does support the use of a PBI threshold similar to the one we identified; the analysis of soil test calibration data for Colwell P showed that calibration curves above and below a threshold PBI of 55 differed (Neuhaus et al. 2018).

The calibration curve we identified for DGT-P is close to those derived in previous work. The DGT-P value at 90% RY for wheat trials conducted between 1968 and 2008 was 34 μg L⁻¹ (n = 131) for soils other than Calcarosols (Speirs et al. 2013). In the high-rainfall zone of south-east Australia, the DGT-P value at 90% RY for wheat in trials 2015–2018 was 25 μg L⁻¹ (n = 12) (McCaskill et al. 2020). It was 41 μg L⁻¹ in the current study but was derived from 10 sites only (out of 40 included in this study). The critical level for DGT-P from the current study was derived from a data set only slightly above the minimum recommended (n = 8) to derive a calibration curve (Dyson and Conyers 2013). A collation and analysis of a larger set of P response data from current cropping systems from sites with a range of soil P buffering levels, paired with measurements of soil P forms, is required to develop DSS capacity using DGT-P.

Objective 2. Investigation of how climatic conditions influence relative yield

The influence of rainfall prior to sowing on RY detected in this study has been observed previously. In a time of sowing x P rate study in South Australia (Mason and McDonald 2021), RY showed a similar relationship with April rainfall for the first time of sowing (26 April−2 May) to that identified here; there was a threshold of about 15 mm after which RY was higher. However, this relationship was not apparent for the second or third time of sowing, indicating that time of sowing influences crop response to P fertiliser via the relative timing of soil P supply and crop P demand (Conyers et al. 2020; Gherardi et al. 2022). Testing for an interaction between time of sowing and rainfall distribution on RY could not be done in this study since there was only one time of sowing at each site. A better understanding of the interaction between time of sowing (a key management decision for WA grain producers (Fletcher et al. 2016)) and rainfall prior to sowing on relative yield for P is needed to inform tactical P fertiliser decisions at sowing.

The impact of rainfall prior to sowing on relative yield is attributed to a combination of factors. The threshold of 8 mm rainfall in a 4-week period may be linked to a minimum water content for mineralisation of organic phosphorus, or improved germination and early growth. Soil carbon mineralisation was detected at gravimetric soil water content as low as 7% in a silt loam after a 1-week incubation (Qiu et al. 2022); however, this threshold has not been examined in sandy soils. Rainfall prior to sowing may also have an additive effect with rainfall around sowing and crop emergence, leading to higher soil water content and faster emergence (e.g. Rinaldi et al. 2005), better early vigour and greater uptake of soil P (McBeath et al. 2012). However, some caution needs to be applied to the use of the rainfall threshold identified in the regression tree analysis. The mean predictive error for monthly rainfall in the SILO database is in the range of 0–5 mm (Jeffrey et al. 2001), about half of the amount of the identified threshold.

Objective 3. Examine whether there were interactive effects between rainfall and soil properties on yield response to P fertiliser

We did not detect a relationship between Colwell P and RY in the current study, but we did detect an effect of rainfall and soil pH_Ca. Our inference is that for soils with a history of P build-up through crop fertilisation, the factors that influence the availability of the existing soil P have a greater influence on crop access to soil P than available soil P levels per se.

How do the outcomes of this study relate to the current knowledge?

The outcomes of the current study differ from previous work in WA, which is attributed to a build-up of soil P in soils with ‘extremely low’ to ‘very low’ PBI. Most fields used for crop production in WA have had a positive P balance for decades; surveys of farmer practice have shown that the mean P balance for crops was 2 kg P ha⁻¹ year⁻¹ (NLWRA 2001; Harries et al. 2021). Long-term positive P balances have led to a build-up of soil P (Weaver and Wong 2011). Soil test calibration curves for Colwell P for P response trials completed in WA from 1958 to 994 (187 trials, r-value = 0.65) and from 1995 to 2011 (107 trials, r-value = 0.57) (Bell et al. 2013) were independent of PBI values, whereas in the current study, we obtained them only for soils with PBI (0–30 cm) values of 9–17. The difference in the capacity to identify a soil test calibration curve for Colwell P may be related to the range of Colwell P levels measured. In the study by Bell et al. (2013), a significant proportion of samples had a Colwell P (0–10 cm) value values of <10 mg kg⁻¹, whereas in the current study, the lowest Colwell P (0–10 cm) was 14 mg kg⁻¹. It may also be due to the over-estimation of plant available P using the Colwell method, due to the extraction of non-labile P (Mason et al. 2013), which is of greater significance in the current compared to the historical work due to the accumulation of soil P precipitated to Al, Ca or Fe in soils used for cropping (McLaughlin et al. 2011).

Outlook

The influence of rainfall prior to sowing detected in this study and sowing date in recent work on wheat yield response to P requires a reassessment of how DSS are structured. At present, DSS used for P fertilisation in Australia are structed around the Mitscherlich equation (e.g. Bowden and Bennett 1974), using site details to parameterise the asymptote (rainfall-limited yield) and the yield with nil P applied (rainfall-limited yield × RY). Rainfall in the period before sowing and growing season rainfall varies greatly from year to year in rainfed cropping systems. Future modelling efforts need to be able to account for the variability in rainfall and management factors such as sowing date on yield response to P fertiliser to inform tactical decision making.

Long-term positive P balances and the resulting build-up of soil P introduces another objective into P fertiliser decisions; managing environmental risk. Phosphorus runoff was detected in WA where grain and livestock production were the dominant land uses, with the source of the P runoff attributed to areas where Colwell P was above critical levels (Sharma et al. 2021). Runoff and leaching losses have not been characterised for the grain-producing region in WA, and this knowledge gap needs to be addressed to link to environmental risk indices (e.g. Weaver et al. 2023) that can be incorporated into P fertiliser decisions for crop production.

The presence of ironstone gravels in the landscapes used for grain production in WA (Holmes et al. 2021) has been a source of concern for P management (e.g. Weaver et al. 2022). In some soil types, ironstone gravel can account for 50–60% of soil weight (e.g. Scanlan et al. 2018), potentially reducing the supply of plant-available P. Previous field research in WA has discarded the gravel fraction when preparing soil samples for chemical analysis so it has not been possible to examine how gravel content influences yield response to P fertiliser (Bell et al. 2013; Weaver et al. 2022). In this study, we did not detect a consistent effect of gravel content on crop response to P fertiliser. However, for future work, the inclusion of detailed measurements of gravel content is still justified because of the extent of soils with gravel present; 15% of soils in the 14 M ha used for grain production in WA have a layer of ferric nodules (ironstone gravel) within the top 15 cm of soil (Geographic Information Services 2016; Holmes et al. 2021).

While the current study has provided valuable advances in our understanding of P response in current farming systems, the relatively small size of the data set poses some limitations. Previous studies have utilised data sets of several hundreds to thousands of trial years (Kuchenbuch and Buczko 2011; Bell et al. 2013). It is possible that some important interactions have not been detected in the current study because of insufficient data, or missing measurements, to capture them, for example a wide range of sowing dates; alternatively, the statistical methods used here could not detect subtle interactions. A larger data set would allow more sophisticated statistical techniques to be used to develop models that could form the basis of DSS (Correndo et al. 2021) and allow a more rigorous analysis using the sliding window approach (Van de Pol et al. 2016).