Study area and soil sampling

Eight silage corn and 23 blueberry fields located in south coastal BC were considered for this study. The soils in the area belong to the order Brunisolic, Gleysolic, Luvisolic, Podzolic, and Regosolic (Soil Classification Working Group 1998). The minimum and maximum general soil properties across all sites were soil pH: 4.9–6.0, total carbon (C): 0.94%–11.35%, total N: 0.07%–0.82%, and Mehlich-3-extractable P: 3.42–195.32 mg·kg⁻¹ (Messiga et al. 2020). Phosphorus fertilization in the silage corn fields included applications of 35 kg P·ha⁻¹ as starter and manure P (Messiga et al. 2020). Phosphorus fertilization recommendations in BC blueberry production guides are based on leaf analysis and are well beyond crop removal (BC berries production guide 2018b). For example, P fertilizer recommendations for a mature blueberry planting range from 20 to 30 kg P·ha⁻¹ for a plant that only exports a maximum of 6 kg P·ha⁻¹·yr⁻¹ including fresh fruit and pruning (BC Berry Production Guide 2018b).

Soil samples (0–15 cm) were collected in silage corn fields in mid-September after harvest in 2018 and 2019 and in the blueberry fields between June and July 2019. In the silage corn fields, 10 soil cores were randomly collected per plot (10 m × 5 m). In the blueberries, each field (approximately 2.0 ha) was visually divided into four blocks, the layer of sawdust on top of the soil was discarded, and 10 soil cores were randomly collected per block. All soil cores were composited on-site, then air-dried, and passed through a 2 mm sieve before analysis. In total, 284 composite soil samples were collected (8 silage corn fields × 24 plots and 23 blueberry fields × 4 blocks).

The average daily temperature ranges from 3.2 °C in December to 18.8 °C in August [30 yr (1981–2010) average]. The local climate is moderate oceanic, characterized by warm, rainy winters and relatively cool, dry summers with mean annual rainfall of 1689 mm [30 yr (1981–2010) average], 261.9 mm of which falls between May and July (Environment and Climate Change Canada 2020).

Phosphorus sorption experiments

The composite soil samples taken from the corn fields were used for P sorption experiments (Messiga et al. 2011). Blueberry sites were not considered for the sorption experiment given the high number of sites and the workload associated. In brief, solutions with a range of P concentrations were formulated in a matrix of distilled water and 100 μL of toluene (to inhibit microbial activity) for batch equilibration with the soil. Initial equilibrating solution P concentrations were 0, 5, 10, 20, 30, 40, 50, 100, 200, and 300 mg P·L⁻¹ as KH₂PO₄. Sorption batches were replicated three times and consisted of 2.0 g air-dried soil mixed in 50 mL centrifuge tubes with 30 mL of the equilibrating solution. The tubes were gently shaken for a 24 h contact time end-over-end at room temperature, centrifuged at 3500 r·min⁻¹ for 10 min and then filtered through Whatman Grade 42 filter paper. An aliquot of the supernatant solution was used to determine P in solution. Phosphorus in the extracts was measured by the molybdate blue colorimetric method which involves a redox reaction in which phosphate ions react to form phosphomolybdate before reading at a wavelength of 880 nm (Murphy and Riley 1962).

Soil analyses

Soil pH was measured in distilled water with a 1:2 soil:solution ratio (Hendershot et al. 1993). Total C and N were determined by dry combustion with a LECO CNS-1000 (LECO Corp., St. Joseph, MI, USA). The Pw was determined by shaking 2.0 g of air-dried soil mixed in 50 mL centrifuge tubes with 20 mL distilled water for 1 h, centrifuging the mixture at 3500 r·min⁻¹ for 10 min, and filtering through Whatman Grade 42 filter paper (Self-Davis et al. 2000). Phosphorus in the extracts was measured by the molybdate blue colorimetric method (Murphy and Riley 1962). Acid ammonium oxalate-extractable P (P_Ox), aluminum (Al_Ox), and iron (Fe_Ox) were determined according to Ross and Wang (1993). Mehlich-3-extractable P (P_M3), aluminum (Al_M3), and iron (Fe_M3) were determined by shaking 2.5 g of air-dried soil with a 25 mL of Mehlich-3 solution (pH 2.3) for 5 min (Mehlich 1984). The concentrations of P_Ox, Al_Ox, Fe_Ox, P_M3, Al_M3, and Fe_M3 were assessed with an inductively coupled plasma optical emission spectrometer (ICAP 7000 series, Thermo Scientific).

The DPS was calculated as the molar ratio of acid ammonium oxalate extracts as follows (van der Zee et al. 1987):

(1)

where P_Ox, Al_Ox, and Fe_Ox are quantified in acid ammonium oxalate extracts (mmol·kg⁻¹). The sum of Al_Ox and Fe_Ox is termed P sorption capacity (PSC, mmol·kg⁻¹). The term α_m is the maximum sorption coefficient (van der Zee et al. 1987; Breeuwsma and Silva 1992). We used an average α_m value of 0.5 to calculate DPS (Benjannet et al. 2018).

The PSI was related to Mehlich-3 extracts as follows (Khiari et al. 2000):

(2)

where P_M3 and Al_M3 are quantified in Mehlich-3 extracts (mmol·kg⁻¹).

Statistical analyses

Sorption data for each silage corn site were fitted using Langmuir and Freundlich isotherms. Sorption isotherms were parameterized using PROC NLIN of SAS version 9 (SAS Institute Inc. 2010). The Langmuir isotherms resulted in the best goodness of fit for all sites. The Langmuir isotherm is described as follows:

(3)

where S (mg P·kg⁻¹) is the sorbed P onto the soil after 24 h contact, S_max (mg P·kg⁻¹) is the P sorption maximum at the highest P concentration, k_a (L mg·P⁻¹) is the binding energy of P, and C (mg P·L⁻¹) is the equilibrium P concentration in solution (Messiga et al. 2011). The goodness of fit of the sorption isotherms to experimental data was assessed by the R² and standard error of the mean.

Pearson’s correlation analysis was performed to assess the linear relationships between S_max and selected soil properties using PROC CORR of SAS (SAS Institute Inc. 2014). To determine environmental risk limits, we investigated the relations among DPS, PSI, and Pw. We used the environmental risk limit for DPS = 25% suggested in the Netherlands (van der Zee et al. 1987) to derive the corresponding critical limits for Pw and then PSI. Finally, we plotted PSI vs. Pw values to verify the consistency of the critical risk limits obtained by using DPS. This approach has been used in several studies in North America including Canada (Khiari et al. 2000; Sims et al. 2002; Benjannet et al. 2018). Graphical and numerical approaches using R² and standard error of the means were the basis to assess the regressions.

We further used these regressions [y = f(x), where x = PSI and y = Pw] to simulate the trends of Pw concentrations for a wide range of soils found in the southern coast of BC. To achieve this goal, we considered a different dataset (pH, P_M3, Al_M3, and Pw) retrieved from a soil survey conducted in 2005 across the southern coast of BC (Table 1; Kowalenko et al. 2007). The 2005 soil dataset included five soil groups (Brunisolic, Gleysolic, Luvisolic, Podzolic, and Regosolic) for corn in rotation with grass and three soil groups (Brunisolic, Gleysolic, and Podzolic) for blueberry production systems. In total, 153 fields were considered from the dataset, including 133 for corn in rotation with grass and 20 for blueberry production systems. The P_M3 and Al_M3 from the 2005 soil survey data were determined by the Mehlich-3 method and were used to calculate corresponding observed PSI values (Mehlich 1984). The Pw concentrations from the survey data were determined according to Self-Davis et al. (2000) and were considered as observed Pw concentrations. The analysis of the observed Pw concentrations was performed using an inductively coupled plasma optical emission spectrometer (ICAP 7000 series, Thermo Scientific) instead of colorimetric method. Considering that P concentrations analyzed by inductively coupled plasma (ICP) are generally higher than colorimetric values by 15%–20% (Adesanwo et al. 2013), we considered a scenario with 0.8 × measured Pw concentration. This difference is generally ascribed to the presence of organic P, which is included in ICP determination but not in colorimetric (Adesanwo et al. 2013). Our assumption was that the regressions (PSI vs. Pw) obtained in the present study would fit the observed PSI values and Pw concentrations collected in 2005. In other words, the regressions will simulate observed Pw concentrations using observed PSI values obtained in 2005. We evaluated the robustness of the regressions by plotting observed Pw (y_observed) concentrations of the 2005 dataset against simulated Pw concentrations (y_simulated). The R², the slope and intercept of the linear regressions, and the root means squared deviation (RMSD) were used to assess the match between y_observed and y_simulated values (Piñeiro et al. 2008). The RMSD was calculated as follows:

(4)

Table 1.

General properties of soils^a (0–15 cm) collected during a provincial soil survey conducted in south coastal British Columbia (Canada) in 2005.

A RMSD value of zero corresponds to a perfect match between simulated and observed Pw concentrations. The closer the RMSD value is to zero, the more accurate the model is.

Phosphorus sorption characteristics of the soils

The experimental data of the sorption study for all eight corn sites were best fitted by the Langmuir equation, with R² values ranging between 0.90 and 0.99 (Fig. 1; Table 2). The S_max values varied from 476 mg P·kg⁻¹ under site 7 to 1914 mg P·kg⁻¹ under site 4 (Table 2). The binding energy of P onto sorption sites (k_a value) also varied among sites, from 0.008 L mg·P⁻¹ for site 1 to 0.041 L mg·P⁻¹ for site 5 (Table 2). Pearson’s correlation analyses showed that S_max values were significantly related to total C (0.62; P = 0.001), PSC (0.42; P = 0.015), DPS (−0.43; P = 0.015), and Al_M3 (0.59; P = 0.001) (Table 3) for corn/grass sites.

Fig. 1.

Phosphorus (P) sorption curves for soils collected between 2018 and 2019 in eight silage corn sites across south coastal British Columbia (TC = total carbon; P_M3 = Mehlich-3-extractable phosphorus; Al_M3 = Mehlich-3-extractable aluminum).

Table 2.

Phosphorus (P) sorption characteristics derived from the Langmuir equations of soils collected in eight corn sites (0–15 cm) across south coastal British Columbia (Canada) in 2018 and 2019.

Table 3.

Pearson’s correlation coefficients of S_max and properties of soils collected in eight corn sites (0–15 cm) across south coastal British Columbia (Canada) in 2018 and 2019.

Relationship between DPS and Pw

The DPS varied widely, ranging between 8.2% and 41.5% among the eight corn sites (Fig. 2a) and between 4.2% and 26.8% among the 23 blueberry production systems (Fig. 2b). The Pw varied between 0.0 and 12.8 mg·L⁻¹ (Fig. 2a) under the corn production systems and between 0.0 and 6.2 mg·L⁻¹ (Fig. 2b) among the blueberry sites. The relationship between DPS and Pw was described by quadratic regressions with R² = 0.89 for the corn (Fig. 2a) and R² = 0.86 for the blueberry sites (Fig. 2b). We pooled the two cropping systems, and the relationship between DPS and Pw was also described by a quadratic function with R² of 0.85 (Fig. 2c). The Pearson’s correlations (0.92 across corn and blueberry sites) derived from the R² show a strong correlation between DPS and Pw. By using this single quadratic regression and considering a critical DPS of 25% (van der Zee et al. 1987), above which the risk of P loss by runoff water is high, we calculated an average critical Pw concentration of 3.7 mg·L⁻¹ (Fig. 2c).

Fig. 2.

Relationships between degree of phosphorus (P) saturation and water-extractable P for (a) silage corn sites (median soil pH = 5.5), (b) blueberry sites (median soil pH = 4.7), and (c) all sites in south coastal British Columbia.

Critical PSI estimation

The PSI varied between 2.2% and 28.6% among corn sites and between 0.62% and 26.2% among blueberry sites (Fig. 3). The PSI were related to DPS through a quadratic regression with R² = 0.83 for the corn sites and a power regression with R² = 0.84 for the blueberry sites, indicating strong correlations (Fig. 3). To refine our estimates of critical PSI, we related PSI and Pw across the two cropping systems. The relationship between PSI and Pw was described by a linear regression with R² = 0.95 for corn sites and a quadratic regression with R² = 0.90 for blueberry sites, indicating strong correlations (Fig. 4). An approximation of the linear regression line shows that, for a critical Pw concentration of 3.7 mg·kg⁻¹, the critical PSI value is 10.4% for the corn sites and 18% for the blueberry sites (Fig. 4). Similarly, if we consider a Pw concentration of 2.0 mg·kg⁻¹, which is considered the lower limit of the risk class, the corresponding PSI is 7.1% for the silage corn sites and 14.4% for the blueberry sites (Fig. 4). Furthermore, a Pw concentration of 5.9 mg·kg⁻¹, which is considered the upper limit of the risk class, corresponds to a PSI of 14.4% for the corn sites and 22.0% for the blueberry sites (Fig. 4).

Fig. 3.

Relationships between degree of phosphorus (P) saturation and P saturation index for silage corn (median soil pH = 5.5) and blueberry sites (median soil pH = 5.5) in south coastal British Columbia.

Fig. 4.

Relationships between phosphorus (P) saturation index and water-extractable P for silage corn (median soil pH = 5.5) and blueberry sites (median soil pH = 5.5) in south coastal British Columbia.

Comparison between observed and simulated Pw values

The PSI in the soil survey dataset varied between 0.65% and 39.6% for the corn sites and between 2.48% and 41.70% for the blueberry sites. The Pw concentrations varied between 1.86 and 29.00 mg·L⁻¹ for the corn sites and between 3.41 and 36.9 mg·L⁻¹ for the blueberry sites. The PSI were related to Pw concentrations by a linear regression with R² = 0.72 for the corn sites (Fig. 5a) and a quadratic regression with R² = 0.93 for the blueberry sites, indicating strong correlations (Fig. 5b).

Fig. 5.

Phosphorus (P) saturation index vs. measured or simulated water-extractable P for (a) silage corn sites (median soil pH = 5.5) and (b) blueberry sites (median soil pH = 5.5) in south coastal British Columbia.

Simulated Pw concentrations paralleled observed concentrations, but the latter were higher than the former for the corn (Fig. 5a) and blueberry sites (Fig. 5b). The RMSD values for the different set of comparisons (simulated vs. observed Pw) were 5.9 for the corn sites (Fig. 6a) and 7.4 for the blueberry sites (Fig. 6b). For the scenario with 0.8 × measured Pw concentration, the RMSD value for the comparisons (simulated vs. 0.8 × observed Pw) was 3.7 for the corn and blueberry sites (Figs. 6a, 6b). In addition, the slopes of the regression lines around a 1:1 line were improved with the new simulations (0.8 × observed Pw values), with 0.71 under the corn sites (Fig. 6a) and 1.21 under the blueberry sites (Fig. 6b) compared with 0.60 and 0.61, respectively, with the original data (measured Pw). This second set of comparisons (Figs. 6a, 6b) highlights the differences between colorimetric and ICP in the determination of P in soil extracts and, therefore, suggest the importance for standardized methods.

Fig. 6.

Measured vs. simulated water-extractable phosphorus (P) for (a) silage corn sites (median soil pH = 5.5) and (b) blueberry sites (median soil pH = 5.5) in south coastal British Columbia.

Phosphorus sorption characteristics

Much research on P sorption characteristics of agricultural soils emphasizes the dominant role of Al, Fe, Ca, and organic matter (OM) content and other soil properties, all of which can increase the capacity of soils to fix newly applied P and thus reduce the efficiency of P fertilizers. The present results indicate that, among the variables analyzed, Al and OM content had the greatest influence on P sorption for the studied soils. The S_max values of the eight soils were in the range of those presented in other studies (Messiga et al. 2011; Wang et al. 2012). However, the wide range of S_max observed indicates differences among soils in their capacity to retain or fix newly applied P from fertilizer and manure. The Pearson’s correlation coefficient between the S_max and Al_M3 was 0.46, and the Pearson’s correlation coefficient between the S_max and OM content was 0.59 (Table 3). The close association between S_max and Al is further highlighted by the significant and positive correlation between S_max and PSC (Al_Ox + Fe_Ox) confirming the role of Al and Fe oxyhydroxides on P sorption characteristics (Table 3). Yan et al. (2013) used a path analysis (a form of multiple regression used to evaluate causal models) to partition Pearson’s correlation coefficients into direct and indirect effects, and their results showed that Al_Ox had a significant direct effect on S_max, whereas the direct effects of Fe_Ox, OM content, and pH on S_max were not significant. The authors concluded that Al_Ox was a dominant contributor to S_max compared with other soil properties (Yan et al. 2013). In the present study, a positive correlation between OM content and S_max could be partly explained by the use of solid manure by dairy farmers as source of nutrients. In a previous study assessing the sorption characteristics of aggregate size particles from no-till and conventional till soils, we discussed three mechanisms through which OM influences S_max (Messiga et al. 2011). Some mechanisms decrease soil P sorption capacity due to direct competition between organic molecules and ligand exchange sites that would otherwise precipitate P (Messiga et al. 2011).

Critical agri-environmental P indicators

The PSC provides an approximation of the quantity of sorption sites in soils available to fix newly applied P (Renneson et al. 2015). A soil with saturated PSC easily releases P to plants during the growing season but also retains less P following rainfall or snowmelt events as they occur during the non-growing season in south coastal BC. The significant relationship (P < 0.001) and R² of 0.85 demonstrate that DPS can be used to predict Pw values across silage corn and blueberry production systems in south coastal BC. Similar relations describing how DPS relates to Pw were obtained across a wide range of soil types in other regions (Pellerin et al. 2006; Renneson et al. 2015; Benjannet et al. 2018). The average critical Pw of 3.7 mg·kg⁻¹ (DPS = 25%) obtained in the present study is close to the environmental threshold of 4.0 mg·kg⁻¹ across Canadian farmlands suggested by van Bochove et al. (2012), but it is lower than the 9.7 mg·kg⁻¹ suggested for soils from Quebec (Pellerin et al. 2006) and the Atlantic provinces of Canada (Benjannet et al. 2018). The method described by Self-Davis was used in the present study and that of van Bochove et al. (2012). The method described by Sissingh (1971) was used by Pellerin et al. (2006) and Benjannet et al. (2018), which likely explains the large difference between Pw concentrations. The Self-Davis method uses a water:soil extraction ratio of 10:1 and shaking time of 1 h, whereas the Sissingh method uses a water:soil extraction ratio of 60:1 and a shaking time of 22 h.

The primary focus of this study was to understand how PSI relates to well-established and widespread environmental indicators such as DPS and Pw in acid to slightly-acid soils of south coastal BC. Two regressions with average R² of 0.84 described the relationships between DPS and PSI, one for the silage corn sites and the other for the blueberry sites (Fig. 3). Similarly, two regressions described the relationships between Pw concentrations and PSI, one for the silage corn sites, with R² of 0.95, and the other for the blueberry sites, with R² of 0.90 (Fig. 4). Interestingly, the critical Pw concentrations of 3.7 mg·kg⁻¹ corresponded to critical PSI values of 10.4% for the silage corn sites and 18.0% for the blueberry sites (Fig. 4). These critical PSI values suggest that the risk of P transport by runoff is high at PSI > 10.4% for silage corn fields and PSI > 18.0% for blueberry fields. The silage corn sites differed from the blueberry sites mainly in their soil pH. Descriptive statistics of soil pH values showed that median soil pH was 5.5 (25% percentile =5.08 and 75% percentile = 5.72) across the silage corn sites and 4.7 (25% percentile = 4.42 and 75% percentile = 5.1) across the blueberry sites (Fig. 5). In a recent study, Penn et al. (2018) showed that P_M3 extraction is decreased with increasing soil pH, partially because of changes in P forms, partial neutralization of extractant pH, and consumption of extractant fluoride by non-P-containing calcium minerals. We found that the occurrence of low pH was also associated with increased Al_M3 concentrations (Fig. 7), which is consistent with other studies (Benjannet et al. 2018). It is, therefore, possible that at low soil pH, Mehlich-3 extraction solution targets forms of P associated with non-oxyhydroxide Al, which overestimates the solubility of labile P. This would explain the higher threshold of critical PSI in the blueberry sites compared with the silage corn sites (Fig. 4). The combination of critical PSI and Pw concentrations as agri-environmental indicators represents an improvement in the way the risk of P loss is assessed among production systems in south coastal BC. The PSI alone can overestimate the risk of P loss to the environment, but coupled with soil pH or Pw, it provides an assessment of the environmental risk of the soil. In silage corn sites, application of manure and supplement of starter fertilizer blends as the main source of N has little effects on the acidity of the soil (Zhang et al. 2018; Messiga et al. 2020). Lower soil pH, in the blueberry sites compared with the silage corn sites is due to the use of ammonium sulfate and other acidifying N fertilizers by local farmers as source of N to inhibit the nitrification process because blueberry plants preferentially absorb ammonium over nitrate (Messiga et al. 2018). Low soil pH also prevails in newly established and young blueberry fields because elemental sulfur is applied to the soil at establishment to decrease the pH, as blueberry plants grow best when soil pH is between 4.2 and 5.5 (Poonnachit and Darnell 2004). The majority of soils under blueberry sites used in the present study had pH values below 5.1 (75% percentile), but in line with the literature a cut-off pH of 5.5 could be used to differentiate silage corn and blueberry sites. Additional studies using soils collected from blueberry sites will be needed to assess these relationships over the entire range of pH (4.2 and 5.5) prevailing in this perennial cropping system. The literature shows that soils with a range of pH as observed in the present study can be divided into two classes according to their P solubility when extracted with Mehlich-3 solution: high solubility with soil pH <5.5; low solubility with soil pH >5.5 (Benjannet et al. 2018; Penn et al. 2018).

Fig. 7.

Relationship between soil pH and Mehlich-3-extractable aluminum for silage corn (median soil pH = 5.5) and blueberry sites (median soil pH = 5.5) in south coastal British Columbia.

Validation of critical Pw value

The simulations mimicked the trends of Pw concentrations for all soils across the two production systems (Figs. 5a, 5b). These results further support the recommendation of two regressions to describe how PSI relates to Pw concentrations in silage corn and blueberry production systems as influenced by soil pH (Benjannet et al. 2018; Penn et al. 2018). However, observed Pw concentrations were on average larger than simulated values, and RMSD values of 5.9 for silage corn sites and 7.4 for blueberry sites were large enough to indicate a significant mismatch between simulated and observed Pw concentrations (Figs. 6a, 6b). One possible reason for this mismatch is that Pw concentrations in soil extracts from the 2005 soil survey were analyzed by ICP, whereas Pw concentrations in soil extracts for the current study were analyzed by the colorimetric method (Murphy and Riley 1962). The literature shows that Pw concentrations determined by ICP are significantly higher than those determined by colorimetric, and the differences are due to the presence of organic P (Adesanwo et al. 2013). To improve the simulations, we considered a scenario with a ratio of 0.8 (80%) between Pw measured by ICP and Pw measured by colorimetric. Accordingly, Pw concentrations obtained in the 2005 soil survey were multiplied by 0.8 (0.8 × observed Pw concentrations) and plotted against simulated Pw concentrations. The mismatch between simulated and observed values were lessened by considering a 0.8 factor, and RMSD value was on average 3.7 for the two cropping systems. These simulations could be improved by testing additional scenarios because the slope of the regression was 0.71 for corn sites and 1.21 for blueberry sites (Figs. 6a, 6b). However, we anticipate that other factors could affect the outputs of the simulations. For example, differences in OM content among soils from the 2005 survey indicate that different ratios, one for each soil, should be used to transform Pw data obtained by ICP. Soil texture is another factor that could affect the output of the simulations and therefore should be accounted for in the grouping of soils. For future work, training and validation Pw concentrations should be obtained using standard methods.

Agronomic and environmental implications for soils with high-legacy P

The PSI is used in acid soils in several Canadian provinces (Sims et al. 2002; Pellerin et al. 2006; Wang et al. 2012; Benjannet et al. 2018). This work is, to our knowledge, the first attempt in south coastal BC to relate PSI and Pw concentrations in two major production systems, silage corn (median soil pH >5.5) and blueberry (median soil pH <5.5), with the goal of defining a critical P risk indicator.

Two critical PSI values were obtained from the present study, 10.5% for the silage corn sites and 18% for the blueberry sites. These values correspond to a critical Pw concentrations of 3.7 mg·kg⁻¹ across the two cropping systems. Above these critical PSI and Pw concentrations, offsite transport of P in surface or subsurface runoff will likely increase. For soils with PSI below critical values, annual P inputs should be planned to maintain agronomical soil P levels. For soils with PSI above critical values, annual P inputs should be less than outputs, and reliance for plant P nutrition should be placed on legacy soil P. For silage corn, one strategy will consist of limiting the use of starter fertilizer P; another strategy is to find ways to adequately address dairy manure application, which is the biggest constraint because of high animal stocking densities in south coastal BC (Reid and Schneider 2019). We recently demonstrated that P fertilizer recommendations for silage corn can be refined by decreasing starter P rates to “5.0 to 7.5” kg P·ha⁻¹ without decreasing dry matter yield at harvest (Messiga et al. 2020). For blueberry production systems, there is a need for scientific evidence on optimum P fertilizer rates to refine actual fertilizer P recommendations. Over a 6 yr study conducted in Prince Edward Island, Canada, plant growth and berry yields of wild blueberry were not affected by rate of soil applied P fertilizer even in soil with low soil test P (Sanderson and Eaton 2008). Another study conducted in Saguenay-Lac-Saint-Jean (Quebec, Canada) between 2009 and 2012 showed a marginal positive effect of P fertilization on berry yield of wild blueberry (Lafond 2020). These studies refer to wild and lowbush blueberries, which may not be relevant for highbush blueberries found in BC. Therefore, additional work is needed to refine P fertilizer recommendations for highbush blueberries.