Tillage Management Impacts on Soil Phosphorus Variability under Maize–Soybean Rotation in Eastern Canada

: Conservation tillage, including no-tillage (NT), is being used increasingly with respect to conventional tillage (CT) to mitigate soil erosion, improve water conservation and prevent land degradation. However, NT increases soil phosphorus (P) stratiﬁcation, causing P runoff and eutrophication. For sustainable P management, fertilization must be balanced between P sources and actual crop demand. To reduce P losses to the environment, it is important to better understand P spatial variability in NT ﬁelds. Little is known about tillage impacts on ﬁeld-scale P spatial variabi-lity in precision agriculture. This study examines tillage impacts on spatial variability of soil-avai-lable P in a maize–soybean rotation, in two commercial ﬁelds, denoted CT (10.8 ha) and NT (9.5 ha), with the aim of improving P fertilizer recommendations in Eastern Canada. NPK fertilizers were applied to the soils (Humic Gleysols) following local recommendations. Soil samples were collected in fall 2014 in regular 35 m by 35 m grids, at 0–5 and 5–20 cm depths, providing 141 and 134 geore-ferenced points for CT and NT ﬁelds, respectively. Available P and other elements were analyzed by Mehlich-3 extraction (M3), and the P saturation index (P/Al) M3 was calculated. Variability of soil-available P in both ﬁelds ranged from moderate to very high (32% to 60%). A mean (P/Al) M3 of 3% was found in both layers under CT, compared to 8% in the 0–5 cm layer and 6% in the 5–20 cm layer under NT. Relationships between P indices and other elements differed between tillage practices. This study highlights the need to improve P fertilizer recommendations in Eastern Canada.


Introduction
Maize (Zea mays L.) and soybean (Glycine max L.) are dominant crops in Eastern Canada, representing 75% of cultivated land and 90% of agricultural production [1]. From 2014 to 2018, 5.5 million tons of these crops were produced, generating total revenue of US $2 billion for the province of Quebec [1]. Owing to their high economic value, these crops are cultivated by local farmers using a variety of tillage methods, including conventional tillage (CT) and conservation tillage methods such as reduced tillage and no-tillage (NT).
CT, which involves preparing the soil with a moldboard plow or a disc plow [2], improves maize productivity by enhancing its root biomass [3] and grain yield [4]. However, this tillage method deteriorates the soil structure and increases the risk of soil erosion, which can lead to higher soil nutrient losses in crops fields [5,6]. To address these problems, conservation tillage practices including NT have been introduced with the broader aim of promoting economically sustainable farming systems [7,8].
Conservation tillage methods, characterized by minimal soil mixing and disturbance, have been adopted increasingly in recent years, owing to their environmental benefits, such as reduced water and wind erosion [9]. The NT method improves soil physical pro-perties by reducing compaction, soil erosion, runoff and drainage risk [9][10][11][12][13][14]. In addition, NT is recognized as the best agronomical practice for maintaining soil integrity in intensively cultivated and large agricultural fields worldwide [11,14,15]. For instance, it is estimated that approximately 9% (125 million ha) of global arable land was under NT ma-nagement in 2011 [14,16]. Canadian farmland under NT management covered an area of 19.5 million ha in 2016 [17]. In Quebec, farmland area under NT has doubled to more than 69% [18], largely due to local government subsidies paid to farmers for using NT practices between 2009 and 2013. Further increases in the use of NT have been observed in recent years [18], owing to the benefits from the standpoint of agronomic [19] and economic performances [20].
Agri-environmental concerns have nonetheless raised general concerns about notillage agronomic practices, in relation to P stratification [21][22][23], transport and runoff to surface waters [24,25], which may cause P eutrophication in rivers [26]. Sustainable management of P involves balancing fertilizer application between P supply (sources) and actual crop demand for soil P [27]. To achieve this objective, it is important to increase our knowledge of the spatial variability of soil P under NT management for developing more rational and cost-effective P fertilization programs that will enhance yields and reduce P losses to the environment. The spatial variability of soil P can be controlled through two main agronomic strategies: (1) variable rate applications (VRA) and (2) the use of management zones (MZs) in precision agriculture.
Many studies [2,[28][29][30][31][32] have examined the spatial variability of soil P under various tillage systems at different sampling, plot or field scales. Cambouris et al. [33] characte-rized the distribution of Mehlich-3 P (P M3 ) concentrations and agri-environmental P saturation (P/Al) M3 at decimetric scale in NT and moldboard plow (MP) plots in a long-term maizesoybean rotation (>20 yr). They observed high CVs associated with P M3 data in both MP (77% and 63% at 0-5 and 5-20 cm, respectively) and NT plots (46% and 66% at 0-5 and 5-20 cm, respectively). Nevertheless, these authors found that the 2D geospatial model related to tillage was not detected by the nested sampling grid used at plot-scale in their study. They suggested the use of an appropriate geostatistical sampling strategy at field scale. Sun et al. [26] reported that studies on the spatial variability of P at plot or experimental scales allow limited interpretation over a larger field area.
In Eastern Canada, local P recommendations have been developed for uniform P a-gronomic applications [34], based mainly on both soil P availability [(P M3 ) and (P/Al) M3 ] index, regardless of field size, soil P variability or their crop productivity potential. Soil variability and associated potential yield differences in a given field are responsible for the development of large field-scale heterogeneity in the distribution pattern of soil nutrients including P [35,36]. P fertilization based on uniform recommendations may result in over-fertilization in large field areas and also create under-fertilized areas [37]. In Quebec, the critical environmental threshold for soil P corresponds to a (P/Al) M3 value of 8% for fine-textured soils [38,39]. Consequently, a better understanding of the spatial variability of soil P in large fields will help to establish accurate P recommendations and ensure that P is applied at the right rate in the right locations.
To our knowledge, no studies have investigated the impacts of contrasting soil tillage systems (CT vs. NT) on the spatial variability of soil P in commercial fields in Eastern Canada under maize-soybean rotation, the dominant legume-cereal association in North America. The main goal of this study was to evaluate the spatial variability of soil-available P under two contrasted tillage systems for developing accurate future P fertilizer recommendations in Eastern Canada. The specific objectives were: to (1) investigate fieldscale spatial variability of soil P indices-referred to as P M3 and (P/Al) M3 -and other selected soil chemical properties (soil pH H2O , total carbon, Al M3 , Ca M3 , Fe M3 ) by using descriptive statistics and geostatistical tools; to (2) study the relationships between soil P indices and these soil chemical properties based on Spearman correlations and multiple regression equations; and to (3) evaluate the current P fertilizer recommendations using prescription maps of kriged values of (P/Al) M3 . The knowledge gained from this research on the effects of tillage management on the spatial variability of soil P indices will provide significant information for optimizing P fertilization in sustainable agriculture systems.

Site Description and Soil Sampling
The study was conducted near Montreal, in the Montérégie region of the province of Quebec, Canada. The Montérégie region is located in the physiographic region of the St. Lawrence Lowlands, which is characterized by very flat landforms (0-5%) [40]. Two commercial fields were selected to represent two contrasting soil tillage practices. One field, denoted the CT field, located at St. Marc-sur-Richelieu, was managed under CT (45 •  In the CT field, tillage has consisted of one plowing operation (20 cm deep) since 2009 using a moldboard plow in the fall after harvest, followed by disking and harrowing (10 cm deep) each spring before seeding. In the NT field, flat direct seeding was performed from 1994 to 2014 onward, with crop residues (30% of residues) left on the ground after harvest. NPK fertilizer sources were as follows: N fertilizer provided from ammonium nitrate (34-0-0); P fertilizer provided from triple superphosphate (0-46-0); and K fertilizer provided from muriate of potash (0-0-60). For each crop, application rates were applied based on previous soil chemical analysis and following province of Quebec recommendations [34]. Both fields were seeded with maize in May 2014 and harvested in late October 2014. At seeding, the NPK fertilizers were band-applied similarly in both fields under the two tillage systems. Using commercial fertilizers, P and K were both applied at the see-ding stage, while N was split at the seeding stage, and at the 6-8 leaf stage, according to local fertilizer recommendations [34]. For maize, the seeding rows were spaced 75 cm apart in both fields. For soybean, the seeding rows were spaced 30 cm apart in both fields. Crop management, maize and soybean fertilization practices were based on local recommendations for the province of Quebec [34].
Both soils were classified as Orthic Humic Gleysols and presented similar soil series (St-Urbain and Kierkoski) [41]. Soil textures for both fields varied from clay to clay loam. They were classified to the same soil texture group (i.e., fine-textured soils with >30% clay) used in the province of Quebec recommendations [34]. Both fields were poorly drained. Meteorological data from the study region were summarized (Table 1). The mean annual temperature is 6 • C and total precipitation is 973 mm at St. Marc-sur-Richelieu. The mean annual temperature is 5.9 • C and total precipitation is 984 mm at La Presentation. An intensive soil sampling was performed on November 2014 using a grid design with a sampling interval of 35 m × 35 m, providing 141 and 134 georeferenced sampling points for the CT and NT fields, respectively (Figure 1). A composite soil sample made up of four cores was taken within a 1 m radius of each sampling point at two soil depths (0-5 cm and 5-20 cm) using a 0.05 m diameter Dutch auger. A total of 282 and 268 soil samples were collected from the CT and NT fields, respectively. At both locations, soil samples were collected after the corn harvest.

FOR PEER REVIEW 4 of 21
Average --14.9 13.9 An intensive soil sampling was performed on November 2014 using a grid design with a sampling interval of 35 m × 35 m, providing 141 and 134 georeferenced sampling points for the CT and NT fields, respectively (Figure 1). A composite soil sample made up of four cores was taken within a 1 m radius of each sampling point at two soil depths (0-5 cm and 5-20 cm) using a 0.05 m diameter Dutch auger. A total of 282 and 268 soil samples were collected from the CT and NT fields, respectively. At both locations, soil samples were collected after the corn harvest. Soil samples were air-dried, ground and sieved through a 2 mm sieve for soil characterization. Soil pHH2O (1:1 water) was measured according to Hendershot et al. [42]. The soils were extracted using Mehlich-3 solution at a 1:10 soil: solution ratio [43], and the concentrations (mg kg −1 ) of available phosphorus (PM3), iron (FeM3), calcium (CaM3) and aluminum (AlM3) were determined by inductively coupled plasma optical emission spectroscopy (ICP-OES; Model, 4300DV, PerkinElmer, Inc., Shelton, CT, USA). The soil P agrienvironmental indicator (P/Al)M3, which is a ratio (%), was calculated from the soil PM3 and AlM3 concentrations. Total carbon (TC, g kg −1 ) content was determined using an Elementar Vario MAX CN analyzer (Elementar 137 Analysensysteme GmbH, Hanau, Germany).
Two topographic parameters were produced for each field using the digital terrain model derived from light detection and ranging (LiDAR) imagery of the province of Quebec. These parameters are elevation and the topographic wetness index (TWI, where TWI = ln(α/tan(β)) is as described by Beven and Kirkby [44] and Sörensen et al. [45]. The value of TWI depends on α, which represents the upslope area per unit contour length, and β, where tan(β) is the local slope. The spatial resolution of these rasters was 1 m. The CT and NT fields were extracted from the respective LiDAR sheets 31H11NE and 31H11SE, using the "extract by mask" ArcGIS tool (ArcGIS Software version 10.3 [ESRI, Redlands, CA, USA]). Mean elevation and TWI values were generated from each sampling point using the buffer method within a 1 m radius and the "Zonal Statistic as Table"   Soil samples were air-dried, ground and sieved through a 2 mm sieve for soil characterization. Soil pH H2O (1:1 water) was measured according to Hendershot et al. [42]. The soils were extracted using Mehlich-3 solution at a 1:10 soil: solution ratio [43], and the concentrations (mg kg −1 ) of available phosphorus (P M3 ), iron (Fe M3 ), calcium (Ca M3 ) and aluminum (Al M3 ) were determined by inductively coupled plasma optical emission spectroscopy (ICP-OES; Model, 4300DV, PerkinElmer, Inc., Shelton, CT, USA). The soil P agri-environmental indicator (P/Al) M3 , which is a ratio (%), was calculated from the soil P M3 and Al M3 concentrations. Total carbon (TC, g kg −1 ) content was determined using an Elementar Vario MAX CN analyzer (Elementar 137 Analysensysteme GmbH, Hanau, Germany).
Two topographic parameters were produced for each field using the digital terrain model derived from light detection and ranging (LiDAR) imagery of the province of Quebec. These parameters are elevation and the topographic wetness index (TWI, where TWI = ln(α/tan(β)) is as described by Beven and Kirkby [44] and Sörensen et al. [45]. The value of TWI depends on α, which represents the upslope area per unit contour length, and β, where tan(β) is the local slope. The spatial resolution of these rasters was 1 m. The CT and NT fields were extracted from the respective LiDAR sheets 31H11NE and 31H11SE, using the "extract by mask" ArcGIS tool (ArcGIS Software version 10.3 [ESRI, Redlands, CA, USA]). Mean elevation and TWI values were generated from each sampling point using the buffer method within a 1 m radius and the "Zonal Statistic as Table" ArcGIS tool (ArcGIS Software version 10.3 [ESRI, Redlands, CA, USA]).

Statistical and Geostatistical Analyses
Descriptive statistics were analyzed using SAS software version 9.4 (SAS Institute Inc., Cary, NC, USA) [46] to investigate soil pH H2O , P M3 , Al M3 , (P/Al) M3 , Fe M3 , Ca M3 , TC, elevation and TWI and their intensity of variability. The minimum and maximum values, the mean, CV and standard deviation of the mean (SDM) were determined for the 0-5 cm and 5-20 cm soil layers from each field. Data were summarized using SigmaPlot software version 14.5 (Systat Software Inc., Palo Alto, CA, USA). The intensity of the variability of the different soil properties was classified based on the approach of Nolin and Caillier [47] using five classes based on the CV: (1) low (CV < 15%); (2) moderate (15% < CV < 35%); (3) high (35% < CV < 50%); (4) very high (50% < CV <100%); and (5) extremely high (CV > 100%). Spearman correlation coefficients (r) between soil P indices [P M3 , (P/Al) M3 ] and other soil chemical properties (Al M3 , Ca M3 , Fe M3 , TC and soil pH H2O ) were determined for each soil layer using SAS software [46]. Spearman coefficients between soil P indices and topography properties (elevation, TWI) were only determined for each topsoil layer. Data were analyzed using the CORR and Univariate procedures.
To determine the relation between the response variables [(P/Al) M3 and P M3 ] and the explanatory variables (Al M3 , Ca M3 , Fe M3 , TC and soil pH H2O ), a mixed regression model was used with a spherical spatial correlation structure between observations. Since measurements were taken at different spatial points (X i , Y i ) and at two different depths (0-5 and 5-20 cm) at each point, a different spherical spatial correlation was assumed for the two depths. Lettingĉ n andα 0 denote the nugget and the range estimates, respectively, the correlation between any two observations a distance d apart, ρ(d) is equal to The regression models were fitted to each field separately. For each response variable, the Box-Cox methodology was used to transform the response variable. The explanatory variables were added sequentially to the model, with the interaction depth factor evaluated each time using the Bayesian information criterion (BIC) until a minimum value was reached. Thus, the final model consisted of the most important explanatory variables and their interaction with the depth factor, which allowed an equation model to be derived for each depth. A pseudo R 2 value was calculated to evaluate the performance of the model. All analyses were performed using the MIXED procedure in SAS (SAS Institute Inc., Cary, NC, USA, release 9.4) at the 0.05 level of significance. The normality assumption was validated using the Shapiro-Wilk statistic on the normalized residuals of the model.
Geostatistics provide the basis for interpolating the spatial variability of soil properties [48][49][50]. Geostatistical analyses of P M3 (mg kg −1 ), (P/Al) M3 (%), Al M3 (mg kg −1 ), TC (g kg −1 ), Ca M3 (mg kg −1 ), soil pH H2O and Fe M3 (mg kg −1 ) data were performed for each soil layer using the GS+ (version 9) software (Gamma Design Software, LLC., Plainwell, MI, USA) [51]. The spatial dependence or structure of these soil properties was evaluated using isotropic and anisotropic semivariograms. Semivariogram parameters were generated for each theoretical model (spherical, exponential and linear). The corresponding nugget (C 0 ), partial sill (C), sill (C 0 + C) and range values of the best-fitting theoretical model were calculated. The partial sill ratio [C/(C 0 + C)] was used to determine the spatial dependence of the soil properties for the CT and NT fields. Semivariograms with a partial sill ratio of <25%, 25% to 40%, 40% to 60%, 60% to 75% or >75% were considered to have a low, low moderate, moderate, strong moderate or strong spatial dependence, respectively [52]. The range indicates the maximum distance at which sample points are correlated [53]. Spatial variability maps were generated for each soil layer using the block kriging method, with a block size of 1 m × 1 m, and were evaluated using cross-validation analysis (R 2 CV ) [54]. Prescription maps were produced using the kriged values of soil P indices (P/Al) M3 for maize and soybean fields using ArcGIS software version 10.3 (ESRI, Redlands, CA, USA), in order to identify different P agri-environmental classes for Quebec and the corresponding areas in both fields. Lastly, accurate P fertilizer recommendations (kg P 2 O 5 ) were developed for maize and soybean in these two commercial fields. These results were compared to the local uniform P fertilizer recommendations [based on (P/Al) M3

Descriptive Statistics of Topography Properties
In the CT field, elevation values ranged from 19.2 to 19.9 m, with a mean value of 19.4 m. The TWI values ranged from 6 to 13 with a mean value of 9. The CVs for elevation (0.8%) and TWI (13%) were low. In the NT field, elevation values ranged from 34.2 to 35.0 m with a mean value of 34.8 m. The TWI values ranged from 7 to 14 with a mean value of 10. The CVs for elevation (0.3%) and TWI (11%) were low. Similar TWI values were obtained, and the CVs for topography properties were low.

Descriptive Statistics of Soil Phosphorus Indices and Other Soil Chemical Properties
In the CT field, similar mean values for the soil P indices [P M3 , and (P/Al) M3 ] and other soil chemical properties were obtained for both soil layers (0-5 cm and 5-20 cm), with the exception of Fe M3 and Al M3 (Figure 2). Mean values of Fe M3 and Al M3 were, respectively, 8% and 12% greater in the 5-20 cm soil layer relative to the 0-5 cm layer. The CVs ranged from 3% to 59% for both soil layers ( Figure 2). The highest CVs were obtained for soil P indices. The CVs were classified as follows: low for soil pH H2O , TC, Al M3 , Fe M3 and Ca M3 , and high to very high for both soil P indices.
In the NT field, mean values of soil P indices, TC and Ca M3 were higher in the 0-5 cm layer than in the 5-20 cm layer. Conversely, the mean values of soil pH H2O , Al M3 and Fe M3 ( Figure 2) were higher in the 5-20 cm layer than in the 0-5 cm layer. The CVs ranged from 7% to 43% for both soil layers ( Figure 2). The highest CVs were obtained for soil P indices. The CVs were classified as follows: low for soil pH H2O ; moderate for TC, Fe M3 and Ca M3 ; low to moderate for Al M3 ; and moderate to high for the two soil P indices ( Figure 2). The intensity of variation of the soil P indices was higher in CT than in NT for both soil layers.

Geostatistical Parameters of Soil Phosphorus Indices and Other Soil Chemical Properties
In the CT field, soil P indices and the other soil chemical properties (except for TC and Ca M3 ) were fitted with spherical and linear models ( Table 2). Pure nugget models (PNs) were used for TC and Ca M3 in both soil layers. Spatial dependence ratios of soil P indices and the other chemical properties ranged from 23% to 85% for both soil layers. Spatial dependence ratios were classified as follows: low moderate for soil pH H2O ; moderate to strong for Al M3 ; low moderate for Fe M3 ; and low to low moderate for soil P indices ( Table 2) Table 2. Geostatistical parameters of the soil chemical properties for two soil layers (0-5 cm and 5-20 cm) from the conventional and no-tillage fields.  2 Sill ratio (%) = [C/(C 0 + C)] × 100; this ratio measures spatial dependence or structure according to Whelan and McBratney (2000). 3 Distance at which a semivariance becomes constant. 4 Coefficient of determination of cross-validation.
The spatial dependence of soil P was higher for both NT soil layers compared to CT soil layers. Conversely, R 2 CV values were generally lower for the NT field relative to the CT field. Higher R 2 CV values (R 2 CV > 0.6) indicate a good fit for kriged map reliability. Consequently, the R 2 CV values for the CT field indicate a relative fit for map reliabilities for most soil chemical properties, including soil P indices, while a weaker fit was found for map reliabilities of soil chemical properties in the NT field.

Spatial Distribution of Soil Phosphorus Indices and Other Soil Chemical Properties
In the CT field, visual similarities were observed between P M3 and (P/Al) M3 spatial distribution maps. Spatial distribution values of soil P indices were similar in both soil layers (Figure 3a,b,g,h). Spatial distribution values of TC were higher in the 0-5 cm soil layer compared to the 5-20 cm soil layer (Figure 3d,j). Spatial distribution values of Al M3 and Fe M3 were lower in the 0-5 cm soil layer compared to the 5-20 cm soil layer (Figure 3c,e,i,k).
In the NT field, visual similarities were also observed between P M3 and (P/Al) M3 spatial distribution maps (Figure 4a,b,g,h). However, spatial distribution values of soil P indices, TC and Ca M3 were higher in the 0-5 cm soil layer than in the 5-20 cm soil layer (Figure 4a,b,d,f-h,j,l). Spatial distribution values of these soil chemical properties ranged from 35 to 19 mg kg −1 , 5.0% to 12.5%, 19 to 27 g kg −1 and 1600 to 4540 mg kg −1 compared to 35 to 95 mg kg −1 , 2.0% to 9.5%, 13 to 25 g kg −1 , 1600 to 2860 mg kg −1 in the 0-5 cm and 5-20 cm NT soil layers, respectively. The spatial distribution value of Al M3 was similar in both soil layers (Figure 4c,i). The spatial distribution value of Fe M3 was lower in the 0-5 cm soil layer than in the 5-20 cm soil layer (Figure 4e,k).    Spatial distribution maps also revealed visual associations between soil P indices and other soil elements. In the CT field, high soil P indices values correspond to the high TC and Ca M3 concentrations in both soil layers (Figure 3a,b vs. Figure 3d,f; and Figure 3g,h vs. Figure 3j,l). In the NT field, the soil P accumulation regions in the soil layer (0-5 cm) correspond to high Ca M3 concentrations, providing evidence of P M3 -Ca M3 accumulation (Figure 4a,b,f). The highest soil P accumulation regions in the soil layer (0-5 cm) correspond to the lowest Fe M3 concentrations.

Relationships between Soil Phosphorus Indices and Other Soil Chemical Properties
In the CT field, the Spearman correlation coefficients with the highest significant relationships were observed between soil P indices and soil TC, Fe M3 and Ca M3 in both soil layers (Table 3). In the NT field, the Spearman correlation coefficients with the highest significant relationships were observed between soil P indices and soil TC, pH water and Ca M3 . In addition, a significant relationship was observed between soil P indices and Fe M3 in the 5-20 cm soil layer in the NT field (Table 3). Table 3. Spearman correlation coefficients of the soil-available phosphorus indices (P M3 ) and (P/Al) M3 in relation to soil chemical and topography properties for two soil layers (0-5 cm and 5-20 cm) in the conventional and no-tillage fields.
Spatial multiple regression equations were generated for the soil P indices in order to investigate the influence of selected soil chemical properties on the variability of soil P indices (Table 4). In the CT field, 45% to 72% of the variability of soil P indices was explained significantly mainly by TC, Fe M3 , Al M3 and Ca M3 . In the NT field, 28% to 39% of variability of soil P indices was only explained significantly by Ca M3 and Fe M3 .

Phosphorus Fertilizer Recommendations Based on Prescription Maps
In the CT field, kriged spatial values for soil (P/Al) M3 ranging from 1% to 4.5% correspond to two Quebec agronomic classes (0% to 2.5%, and 2.6% to 5%) for maize and soybean. According to these kriged values for (P/Al) M3 , the P fertilizer recommendations should be 60 or 80 kg P 2 O 5 ha −1 for maize and 20 or 60 kg P 2 O 5 ha −1 for soybean (Figure 5a,e). The applied areas (A 1 and A 2 ) from these kriged values were 6.1 and 4.7 ha. Thus, in the CT field, accurate total P recommendations were 770 and 460 kg P 2 O 5 for maize and soybean, respectively (Table 5). Based on the average value of (P/Al) M3 in the CT field (2.7%), the uniform P agronomic recommendations should be 60 kg P 2 O 5 ha −1 , for a total of 648 kg P 2 O 5 for maize, and 20 kg P 2 O 5 ha −1 , for a total of 216 kg P 2 O 5 for soybean (Table 5).    No-tillage 2.6 20 52 6 0 0 0.9 0 0 52 0 1 A: specific area (ha) from each specific management zone measured using ArcGIS. 2 Rate: P specific recommendation (kg P 2 O 5 ha −1 ) corresponding to each kriged value of (P/Al) M3 from the specific area. 3 Q = A × R: this value represents the accurate total P recommendation (kg P 2 O 5 ) applied in each delimited area (A), based on each kriged value of (P/Al) M3. 4 Total P fertilizer recommendation (kg P 2 O 5 ) or sum of (Q 1 + Q 2 + Q 3 ) applied from all specific areas A 1 , A 2 and A 3. 5 Local P fertilizer recommendation (kg P 2 O 5 ) from the mean value of (P/Al) M3 traditionally applied in the field according to the Guide de Référence en Fertilisation du Québec (CRAAQ, 2010).
In the NT field, kriged spatial values for soil (P/Al) 3 , which ranged from 2% to 12.5%, correspond to three Quebec agronomic classes, i.e., 2.6% to 5%, 5.1% to 10.0% and 10.1% to 15% for maize, and 2.6% to 5%, 5% to 7.5% and 7.6% to 15% for soybean. The P fertilizer recommendations should be 20, 40 or 60 kg P 2 O 5 ha −1 for maize, according to these (P/Al) M3 kriged values (Figure 5b). The applied areas corresponding to these kriged values were 0.7, 8.8, 6.9 and 2.6 ha ( Table 5). For soybeans, the P fertilizer recommendations should be 20 or 0 kg P 2 O 5 ha −1 , according to these (P/Al) M3 kriged values (Figure 5f). The applied areas corresponding to these kriged values were 2.6, 3 and 6.5 ha (Table 5). In the NT field, the accurate total P recommendations were 366 and 0 kg P 2 O 5 for maize and soybean, respectively (Table 5). Based on the average value of (P/Al) M3 in the NT field (7.9%), the uniform P agronomic recommendations should be 40 kg P 2 O 5 ha −1 , for a total of 380 kg P 2 O 5 for maize, and 0 kg P 2 O 5 ha −1 , for a total of 0 kg P 2 O 5 for soybean (Table 5).

Effects of Topography Properties on Phosphorus Transport
Both fields used in this study were flat, as confirmed by previous soil studies conducted in the region [35,40]. The TWI is a factor that promotes P transport in temperate agricultural landscapes [55,56]. Mean TWI values were higher than in Li et al. [57]. A higher TWI is associated with increased surface runoff, which may result in greater export of nutrients [57], such as P. However, the CV values of TWI were lower than in other studies [57,58]. Thus, similar TWI values from the contrasting tillage fields and their res-pective low CVs provide evidence of failure to take account of these topographic properties in relation to P transport and spatial variability. This is confirmed by results showing a lack of significant relationships between these topography properties and soil P indices (Table 3).

Effects of Soil Tillage Practices on Soil Phosphorus Stratification
The mean values of soil P indices were higher in the 0-5 cm soil layer relative to the 5-20 cm layer under the NT system. In contrast, the mean values were similar in both soil layers in the CT field. Our results show that P stratification exists (Figure 4a,b,g,h). Soil P stratification is defined as high soil P concentrations at the soil surface and rapidly decreasing concentrations with depth [59]. Many studies [21,22,[60][61][62][63] have found soil P stratification in agricultural plots managed using the NT system relative to the CT.
In this field-scale study, soil P stratification can be explained by crop residues and the lack of mixing of soil and fertilizers. Previous studies [61,64] reported that P stratification in soils under NT has three causes: (1) annual broadcast or band application of P fertilizer at or near the same sowing row; (2) absence of mixing of P fertilizers with soil and crop residues left on the soil surface after harvest; and (3) immediate leaching of P from crop residues.
The results obtained show that increased soil P accumulation is associated with NT, which points to a higher agri-environmental risk for NT field management in Eastern Canada. The mean (P/Al) M3 index increased (7.9%) in the 0-5 cm NT soil layer relative to the 5-20 cm NT soil layer (6%). These results also have environmental implications for soil P. The mean (P/Al) M3 value for the NT field (7.9%) was close to Quebec's national critical threshold for P saturation for fine-textured soils (8%), indicating that (P/Al) M3 values higher than this can result in water contamination. Consequently, NT affects the sustainable management of soil P, increasing the risk of P pollution in Eastern Canada.

Variability of Soil Phosphorus Indices under the Different Soil Tillage Practices
Among all properties measured in both fields, the highest CVs were obtained for soil P indices. These results are in line with other studies [2,65]. In this field-scale experiment, the CVs of soil P indices ranged from moderate to very high (32% to 60%). Many studies [30,32,33] obtained similar results under different tillage systems at different scales. For instance, Malvezi et al. [2] reported CVs (n = 100; 0-20 cm) for soil-available P ranging from moderate to high (32% to 76%) under two contrasting CT and NT experimental plot-scale systems cultivated for more than three decades.
The CVs of soil P indices were high (32% to 43%) in the NT field. Other studies [32,66,67] reported that high CVs of soil P in NT fields were mainly caused by (1) limited mixing of soil, crop residues and fertilizers, and (2) the high residual P levels in soil surface in NT fields. A high CV of soil P in NT fields has implications for developing future agronomic strategies for sustainable management of P based on delineating MZs or VRA for accurate P recommendations in precision agriculture.
The intensity of variation of soil P indices was lower in both soil layers (0-5 cm and 5-20 cm) under the NT system, compared to CT. These results are similar to the findings of Cambouris et al. [33], in Eastern Canada, and stand in contrast to those of Malvezi et al. [2] in southern Brazil. Consequently, variability of soil P may be reduced under NT at the scale used in the present study. The lower CV obtained with the NT system will positively impact future soil P sampling strategies, based on a limited number of soil samples.

Geostatistics of Soil Phosphorus Indices under Different Soil Tillage Practices
The spatial dependence of soil P indices was moderate (23% to 46%), indicating that soil P variability in the contrasting tillage fields resulted from interaction between intrinsic (soil types) and extrinsic soil factors, such as soil/crop management practices (tillage mode and application of fertilizers or manure). Moderate spatial dependence (25% < R < 75%) results from interaction between intrinsic and extrinsic soil factors [48,68]. These results are similar to those of Dalchiavon et al. [30], who found a moderate spatial dependence (54%) for soil P in a 5 ha soybean field under a NT system. In contrast, other studies [2,33] reported pure nuggets (PNs) under CT and NT long-term experiments. These PNs were explained by (1) the plot scale of these studies, and (2) a strong overall variability of soil P indices, which overshadowed soil P patterns at the plot-scale used in the study. Spatial ranges revealed that the sampling grid used to measure spatial variability in both fields was appropriate. The present field-scale study also revealed that the 2D geospatial model related to tillage was detected for soil P, based on the sampling grid used initially, compared to the plot-scale experiment conducted by Cambouris et al. [33]. In this earlier study, the 2D geospatial model related to tillage was not detected by the sampling grid used. Furthermore, the smallest range values, 46 m and 55 m, corresponded to soil P indices from CT and NT fields, respectively ( Table 2), indicating that these minimum ranges will strongly influence the optimum grid for determining soil P and other chemical properties in both fields [69]. Consequently, future strategies for sampling soil P and other chemical properties for geostatistical research should be planned according to the tillage systems used, taking into account the range of P values corresponding to each tillage system. Spatial dependence of soil P was greater under the NT system compared to CT. This can be explained by the reduction in soil tillage activities under NT fields over a long period, which affects the long-term spatial structure of soil P at field-scale. In contrast, the lowest R 2 CV values for soil P indices were observed with stronger spatial structures in the NT field. Based on the hypothesis of Kravchenko [70], we can infer that the geostatistical models traditionally used are not suited to NT fields. Previously, Mallarino [32] confirmed this idea in his field-scale study on NT systems in Iowa. Further studies need to be carried out using other adapted models to achieve better prediction of spatial structure of soil P indices in NT fields.

Spatial Distribution Maps and Relationships between Soil Phosphorus Indices and Other Soil Chemical Properties
The R 2 CV is an important reliability factor for producing kriged spatial maps of soil P indices with the aim of preventing P pollution. R 2 CV values (R 2 CV < 0.2) were obtained for most soil properties. These values are similar to those of other studies [33,36] and lower than those reported by Perron et al. [71]. The lowest R 2 CV values are due to the weak spatial dependence of the soil properties [70]. Furthermore, kriged spatial maps of soil P indices also revealed soil tillage impacts on the sustainable agri-environmental management of P, particularly for identifying zones with agronomic constraints on P fertility, which are associated with the potential loss of P to the environment. For instance, low P fertility values (1% to 4.5%, Figure 3b) were observed in the CT field. High P fertility values (5% to 12.5%) were observed in the NT field (Figure 4b). At the environmental scale, several NT zones had high P values (8% to 12.5%) that exceeded the national threshold of P saturation (8%) for fine-textured soils [38,39]. Consequently, using the NT system may increase the P fertility level and the risk of P pollution. In addition, visual associations from spatial maps of soil P indices and the other elements (TC, Ca M3 and Fe M3 ) confirmed the influence of these different chemical properties on soil P under the two contrasting tillage practices. Recently, Nze Memiaghe et al. [36] reported similar observations for two contrasting grassland systems in Eastern Canada.
As shown in Table 3, the NT system may cause P accumulation in soil and induce significant changes in soil chemical properties. These significant relationships between soil P indices and other chemical properties depend on tillage effects within fields, owing to a decrease in soil pH. In this study, mean pH values from the CT and NT fields were 7 and 6.8, respectively. Piegholdt et al. [64] found similar mean pH values for CT and NT plot soil layers (0-5 cm and 5-25 cm) on Luvisols and Phaeozems in a study on long-term (7 yr) tillage effects. A decrease in soil pH H2O is commonly observed in intensive NT soils compared to CT soils, as reported by Malvezi et al. [2]. This is probably due to the acidifying effects of N fertilizers associated with the decomposition of crop residues [8,[72][73][74][75].
Furthermore, significant relationships were observed between soil pH H2O and both soil P indices (r = 0.23 for P M3 , p < 0.01; r = 0.26 for (P/Al) M3 , p < 0.01), particularly in the 0-5 cm NT soil layer. Similar results were previously reported in gleysols cultivated intensively from this same region [35] and other soils under intensive NT system [15,30].
Other relationships between soil P indices with soil TC, Ca M3 and Fe M3 were observed under both NT soil layers. This was reported by previous studies [15,30,35,64]. Nevertheless, this soil P accumulation in relationship with these soil elements depends also on significant relationships between soil pH H2O and these chemical properties. Thus, significant relationships between soil pH H2O and Ca M3 were observed in the 0-5 cm (r = 0.76; p < 0.001) and 5-20 cm (r = 0.67; p < 0.001) NT soil layers. Significant relationships between soil pH H2O and Fe M3 in these same layers were (r = −0.68; p < 0.001) and (r = −0.73; p < 0.001), respectively. However, no significant relationships between soil pH H2O and TC were observed in the 0-5 cm (r = 0.07; ns) and 5-20 cm (r = −0.1; ns) NT soil layers. Finally, the results from Table 4 reveal that soil tillage practices induced significant changes in TC, Al M3 , Fe M3 and Ca M3 , affecting the variability of soil P indices. Thus, it is necessary to consider the variability of these chemical properties and their respective spatial structures in order to understand soil P indices variability within fields and according to tillage effects, because these soil chemical properties impact significantly on soil P accumulation under NT systems.

Phosphorus Recommendations and Environmental Implications
As observed in the kriged spatial maps, the soil P accumulation in the first 5 cm of the NT field highlights the environmental risk of P pollution. Average (P/Al) M3 values were 2.7% and 7.9% in the CT and NT fields, respectively, after 20 years. Based on current local recommendations [34], P fertilizer applications are 60 kg and 20 kg P 2 O 5 ha −1 for maize and soybean, respectively, in the CT field. In comparison, the P fertilizer recommendations are 40 kg and 0 kg P 2 0 5 ha −1 for these respective crops in the NT field. Nevertheless, these Quebec local P fertilizer recommendations were developed for uniform P agronomic applications, regardless of field size, soil P variability or their crop producti-vity potential. Soil variability and potential yield differences in a given field are responsible for large field-scale heterogeneity in the distribution pattern of soil nutrients, such as P [35,36].
Variability of soil P indices in both fields ranged from moderate to very high (32% to 60%), indicating that the existing uniform P fertilizer recommendations for maize and soybean are not suited to large (10 ha) fields with contrasting tillage. P fertilization based on uniform recommendations may result both in over and under-fertilization in a large field area [37]. Consequently, accurate P fertilizer recommendations were developed for these fields planted with maize and soybean crops, taking into account (1) spatial distribution values for soil P resulting from high soil P variability in these large fields and (2) the corresponding Quebec agronomic P recommendations for these respective crops ( Figure 5).
For instance, accurate P fertilizer recommendations were developed for maize (770 kg P 2 O 5 ) and soybean (460 kg P 2 O 5 ) in the CT field (Table 5). In comparison with local P uniform recommendations (648 kg and 216 kg P 2 O 5 for maize and soybean, respectively), this represents an increase in P 2 O 5 of up to 19% and 113% for maize and soybean, respectively. The P fertilizer recommendation for the NT field consisted of a lower amount for maize (366 kg P 2 O 5 ha) and no application (0 kg P 2 0 5 ) for soybean (Table 5). This represents a reduction of 14 kg P 2 O 5 for maize and none (0 kg P 2 O 5 ) for soybean, compared to the local uniform P recommendations (380 kg P 2 O 5 and 0 kg P 2 O 5 ). As this example shows, a high CV of soil P in both fields underscores the importance of applying P fertilizer at the right rate and in the right location to avoid creating over-and under-fertilized areas. Thus, there is a potential for controlling spatial variability of soil P in order to increase farm field productivity and yield while reducing P environmental losses.

Conclusions
This study aimed to examine the impacts of tillage (conventional and no-tillage) on the spatial variability of soil-available P at field scale in a maize-soybean rotation in Eastern Canada with the ultimate goal of improving P fertilizer recommendations. The results revealed that the soil P indices were similar in both layers (0-5 cm and 5-20 cm) under CT. However, in the NT field, these P parameters were higher in the 0-5 cm layer relative to the 5-20 cm layer. These findings have implications for the sustainable management of soil P. The NT system increased soil P [(P/Al) M3 ] accumulation (7.9%) compared to the CT method (2.7%), presenting a greater risk of P pollution, as observed in the kriged spatial maps. Relationships between soil-available P indices and other chemical properties differed between the contrasting tillage practices.
The variability of soil P indices in both fields ranged from moderate to very high (32% to 60%), indicating that the uniform P fertilizer recommendations currently used for maize-soybean rotations are not suitable for large (10 ha) fields. The use of kriged maps demonstrated the importance of applying P at the right rate in the right location and using this information to develop accurate future P recommendations. The geostatistical models traditionally used to derive soil P indices are not suited to NT fields. Further research should be conducted using other, improved models for better prediction of the spatial structure of soil P indices in NT fields.