Evaluation of Nitrogen Fertilization Systems Based on the in-Season Variability in the Nitrogenous Growth Factor and Soil Fertility Factors—A Case of Winter Oilseed Rape ( Brassica napus L.)

: Application of nitrogen (N) in contrastive chemical form changes availability of soil nutrients, a ﬀ ecting crop response. This hypothesis was evaluated based on ﬁeld experiments conducted in 2015 / 16 and 2016 / 2017. The experiment consisted of three nitrogen fertilization systems: mineral-ammonium nitrate (AN) (M-NFS), organic-digestate (O-NFS), 2 / 3 digestate + 1 / 3 AN (OM-NFS), and N rates: 0, 80, 120, 160; 240 kg ha − 1 . The content of nitrogen nitrate (N-NO 3 ) and available phosphorus (P), potassium (K), magnesium (Mg) and calcium (Ca) were determined at rosette, onset of ﬂowering, and maturity of winter oilseed rape (WOSR) growth from three soil layers: 0.0–0.3, 0.3–0.6, 0.6–0.9 m. The optimum N rates were: 139, 171 and 210 kg ha − 1 for the maximum yield of 3.616, 3.887, 4.195 t ha − 1 , for M-NFS, O-NFS, OM-NFS. The N-NO 3 content at rosette of 150 kg ha − 1 and its decrease to 48 kg ha − 1 at the onset of ﬂowering was the prerequisite of high yield. The key factor limiting yield in the M-NFS was the shortage of Ca, Mg, O-NFS—shortage of N-NO 3 . Plants in the OM-NFS were well-balanced due to a positive impact of the subsoil Mg and Ca on the N-NO 3 content and productivity. The rosette stage was revealed as the cardinal for the correction of WOSR N nutritional status.


Introduction
In current agriculture, the amount of food production depends on the use of N fertilizers, both mineral and organic [1]. In modern agriculture, which is dominated by varieties of a huge yield potential, resulting from a high rate of biomass growth, the requirement for N is high [2]. This requirement can be fulfilled by taking into account the required rate of a crop plant growth and the uptake of the adequate amount of N, and synchronization of its supply with a plant requirement. N is taken up by plants in two chemical forms, i.e., as nitrogen nitrate (N-NO 3 , NN), and ammonium (N-NH 4 ) [3,4]. This N form, as opposed to ammonium, does not undergo fixation by soil colloids therefore being easily available to a crop plant within a broad range of soil pH [3]. Nitrate N, in spite of a higher metabolic cost of assimilation, results in the production of carbohydrates, consequently leading to a higher rate of plant growth with respect to ammonium [5,6]. In the light of present knowledge, NN, in fact its soil resources, can be considered as the soil nitrogenous growth factor (NGF). This hypothesis is also supported by thousands of scientific papers in the area of agriculture

Site Description
The primary sources of data are two series of field experiments with winter oilseed rape (Brassica napus L.), which were carried out during the 2015/2016 and 2016/2017 seasons in Baniewice (53 • 05 N; 14 • 36 E), Poland. The field experiment was conducted on soil with texture loamy sand in the top-soil over sandy loam, classified as Albic Luvisol (World Refernce Base -WRB, 2016) [38]. The content of the available nutrients was measured each year just after a fore-crop harvesting from three soil depths of 0.0-0.3, 0.3-0.6, and 0.3-0.9 m. The content of available nutrients was favorable for high-yielding WOSR ( Table 1).
The local climate, classified as intermediate between atlantic and continental, is seasonally variable ( Table 2). Precipitation during the period extending from January to July amounted to 332 mm in 2016 and to 622 mm in 2017 (237 in July). In June, a critical month with respect to pods and seed growth [39], the amount of rainfall was extremely low. Air temperatures were in both years higher in comparison to the long-term average.  [40]; 2 average ± standard deviation; 3 fertility class: VL-very low; L-low; M-medium; H-high; VH-very high [41].

Experimental Design
The field experiment was arranged as a two-factorial design. The first factor was the nitrogen fertilization system (NFS), related to the type of N fertilizer applied: (i) ammonium nitrate, 34-0-0, (acronym M-NFS), (ii) organic N: digestate (O-NFS), (iii) mixed: 2 / 3 of digestate + 1 / 3 of ammonium nitrate (AN) (OM-NFS). The rate of applied N was: 0 (N control), 60, 120, 180, and 240 kg ha −1 . Each year, digestate was applied at the end of November, and AN on the 1st of March of the following year. The applied digestate on average contained: DM: 72, N-NH 4 : 7.2, P 2 O 5 : 5.45, K 2 O: 2.85, MgO: 1.88, S: 1.99 kg m −3 . Winter barley was a fore-crop for WOSR. A hybrid WOSR variety Impression Cl, characterized by a high yielding potential for medium fertile soils, was used as a test crop. It was sown at the rate of 3.0 kg ha −1 at the end of August and harvested at the end of July from an area of 12 m 2 . The seeds were harvested at maturity when the moisture content was 8% dry weight. Table 2. Main characteristics of meteorological conditions during the winter oilseed rape (WOSR) growing season on the background of the long-term averages 1 . VIII  IX  X  XI XII  I  II  III  IV  V  VI

Chemical Analysis and Indices
Soil samples were collected from each field three times a year, i.e., (i) at the onset of the stem elongation phase of WOSR growth (BBCH 30, STME), (ii) at the onset of flowering (BBCH 60, FL), and (iii) at the crop physiological maturity (BBCH 89). Soil samples at each stage of WOSR growth were taken from three soil depths as follows: 0-0.3, 0.3-0.6, and 0.6-0.9 m.
The soil content of NH 4 -N and NO 3 -N was determined in field-fresh (not air dried) soil samples within 24 h after sampling. Twenty-gram soil samples were shaken for 1 h with 100 mL of 0.01 M calcium chloride (CaCl 2 ) solution (soil/solution ratio 1:5; m/v). Concentrations of NO 3 -N were determined by the colorimetric method using flow injection analyses (FIAstar5000, FOSS). The method of NO 3 -N concentration analysis consists of two basic steps: reduction from nitrate to nitrite by using a cadmium column, followed by colorimetric determination of nitrite based on the Griess-Ilosvay reaction with N-(1-Naphtyl) ethylene-diamine dichloride as a diazotizing agent. The measurement was performed at a wavelength of 540 nm. The total content of nitrate nitrogen (NN) was expressed in kg ha −1 . The N min content for a given soil layer was calculated using indices, which were constructed based on a soil textural class and soil bulk density [42].
Soil pH was measured in a 1 M potassium chloride (KCl) solution (soil/solution ratio 1:2.5; m/v). The content of plant available nutrients, including P and K, Mg, and Ca, was determined by the Mehlich 3 method [39]. The P concentration in the extract was determined by the colorimetric method, using the molybdenum blue method. Concentration of K in the collected extracts was determined by flame photometry. Concentration of Mg and Ca were determined by atomic absorption spectrometry (AAS)-flame type.
The partial productivity of NN was calculated based on yield (Y, kg ha −1 ) and the amount of NN present in a soil layer of 0.0-0.9 m at BBCH 30 (kg ha −1 ): PFP N30 = Y/N 30 , kg seeds kg −1 of N (1)

Statistical Analyses
The effect of individual research factors (year, N fertilization system, and N rate) for each soil layer and the interaction between them was assessed by two-way ANOVA. Differences between the mean values were compared using the HSD -honestly significant difference test, according to the Tukey method, where the significance level was assessed at 0.05. STATISTICA 13 ® software was used to conduct all statistical analyses. In the second step of the diagnostic procedure, the relationships between variables representing soil properties were analyzed by principal component analysis (PCA) (StatSoft, Inc., Krakow, Poland 2013). In the third step of the diagnostic procedure, stepwise regression Agronomy 2020, 10, 1701 5 of 27 was applied to define an optimal set of variables for a given N characteristic. In the computational procedure, a consecutive variable was removed from the multiple linear regressions in a step-by-step manner. The best regression model was chosen based on the highest F-value for the model and the significance of all independent variables.

In-Season Variability in the Nitrate Nitrogen (N-NO 3 , NN) Content
A significant impact of the nitrogen fertilization system (NFS) was recorded in 2017 (Table 3). In this particular year, a significantly higher seed yield was recorded in the organic-mineral (OM-NFS) treatment, in which 2 / 3 of the total N rate were applied in the form of digestate and the remaining one ( 1 / 3 ) as ammonium nitrate. In 2016, a year with a water shortage during WOSR flowering, no significant difference between NFSs was observed, but a higher yield was harvested from the OM-NFS (Tables 2 and 3). The effect of the progressively increased N rates was revealed each year. In both years, the effect of the applied N was only significant with respect to the absolute control (AC). In 2016, the relative yield increase on the plot fertilized with N applied at a rate of 60 kg ha −1 , compared to AC, was 63%, whereas in 2015 it doubled, reaching 115%. Yield, taking into account both experimental factors and year, significantly depended on the interaction of N rates (N) and year (Y × N) and on NFS and N rates (NFS × N). Due to its significance for all the other studied N characteristics, the second interaction was considered as the most representative for this study (Table 3). As shown in Figure 1, an N rate of 60 kg ha −1 , irrespective of the NFS, resulted in a significant yield increase. The further yield response to increasing N rates was NFS specific. For the M-NFS, an N rate of 120 and 180 kg ha −1 did not result in any yield increase. An N rate of 240 kg ha −1 resulted in a yield drop as compared to the N rate of 60 kg ha −1 . A similar trend was observed for the O-NFS, but without any drop with respect to the treatment with the N rate of 60 kg ha −1 . A significantly different trend was observed for the OM-NFS, in which yield increased significantly in accordance with the progressively applied N rate. The regression models of WOSR yield response to the N f rate are as follows: Agronomy 2020, 10, x 6 of 29 Yield, taking into account both experimental factors and year, significantly depended on the interaction of N rates (N) and year (Y × N) and on NFS and N rates (NFS × N). Due to its significance for all the other studied N characteristics, the second interaction was considered as the most representative for this study (Table 3). As shown in Figure 1, an N rate of 60 kg ha −1 , irrespective of the NFS, resulted in a significant yield increase. The further yield response to increasing N rates was NFS specific. For the M-NFS, an N rate of 120 and 180 kg ha −1 did not result in any yield increase. An N rate of 240 kg ha −1 resulted in a yield drop as compared to the N rate of 60 kg ha −1 . A similar trend was observed for the O-NFS, but without any drop with respect to the treatment with the N rate of 60 kg ha −1 . A significantly different trend was observed for the OM-NFS, in which yield increased significantly in accordance with the progressively applied N rate. The regression models of WOSR yield response to the Nf rate are as follows: Y-M-NFS = −0.00009N 2 + 0.025N + 1.88 for n = 5, R 2 = 0.93, P ≤ 0.05 (2) Y-O-NFS = −0.00007N 2 + 0.024N + 1.83 for n = 5, R 2 = 0.94, P ≤ 0.05 Y-OM-NFS = −0.00005N 2 + 0.021N + 1.99 for n = 5, R 2 = 0.97, P ≤ 0.05 (4) The key attributes of quadratic models, such as the optimum N rate (Nop) and the respective maximum yield (Ymax), were used to distinguish the tested NFSs. The lowest Nop of 138.9 kg ha  The NN content at the rosette stage of WOSR growth significantly responded to NFSs in 2016 (Table 3). In this particular year, the highest NN content was recorded in the M-NFS, followed by OM-NFS, and O-NFS. In 2017, the NN content in the M-NFS was almost at the same level as in 2016, but no significant differences between NFS treatments were recorded. As shown in Figure 2 The key attributes of quadratic models, such as the optimum N rate (N op ) and the respective maximum yield (Y max ), were used to distinguish the tested NFSs. The lowest N op of 138.9 kg ha −1 was the attribute of the M-NFS. The respective Y max amounted to 3.616 t ha −1 . The parameters of the O-NFS were 171.4 kg ha −1 of applied N, and the Y max of 3.887 t ha −1 . The highest values of both parameters were found for the OM-NFS. A N op of 210 kg ha −1 resulted in a theoretical Y max of 4.195 t ha −1 .
The NN content at the rosette stage of WOSR growth significantly responded to NFSs in 2016 (Table 3). In this particular year, the highest NN content was recorded in the M-NFS, followed by OM-NFS, and O-NFS. In 2017, the NN content in the M-NFS was almost at the same level as in 2016, but no significant differences between NFS treatments were recorded. As shown in Figure 2, the NN content followed different regression models in response to increasing N rates. The cubic regression model fitted the M-NFS and O-NFS data the best. The course of the model for the M-NFS showed a significant response to the N rate of 60 kg ha −1 , next stabilizing up to the N rate of 180 kg ha −1 . A sudden increase in the NN content was recorded following the N rate of 180 kg ha −1 . The course of this model suggests huge NN resources during STME in the plot fertilized with 240 kg ha −1 of N. The course of NN on the O-NFS was completely different. Its content increased exponentially up to an N rate of 145.7 kg ha −1 , then significantly decreased. This course, in contrast to the M-NFS, indicates a full exploitation of NN resources in the plot fertilized with 240 kg ha −1 of N. The NN course for the OM-NFS fitted the quadratic function the best, and was characterized by an exponential increase in the content of NN from the N rate of 6.9 kg ha −1 , thus indicating huge NN resources during the STME in plots fertilized with the highest N rates.
Agronomy 2020, 10, x 7 of 29 sudden increase in the NN content was recorded following the N rate of 180 kg ha −1 . The course of this model suggests huge NN resources during STME in the plot fertilized with 240 kg ha −1 of N. The course of NN on the O-NFS was completely different. Its content increased exponentially up to an N rate of 145.7 kg ha −1 , then significantly decreased. This course, in contrast to the M-NFS, indicates a full exploitation of NN resources in the plot fertilized with 240 kg ha −1 of N. The NN course for the OM-NFS fitted the quadratic function the best, and was characterized by an exponential increase in the content of NN from the N rate of 6.9 kg ha −1 , thus indicating huge NN resources during the STME in plots fertilized with the highest N rates. The NN content at BBCH 60, i.e., at the onset of WOSR flowering (FL), was significantly driven by both experimental factors and the interaction between them. It is necessary to stress that the NN content at end of the STME was significantly correlated with its status at the beginning of this phase, i.e., BBCH 30 (Table A1). In 2016, the impact of NFS on the NN content, evaluated independently of N rates, was the same as at BBCH 30. In 2017, a significantly higher amount of NN was recorded in the M-NFS. The effect of the progressively increased N rates was significantly different in both years. In 2016, the NN content increased up to the N rate of 120 kg ha −1 and then fell significantly, whereas in 2017 it increased in accordance with the increased N rates. As shown in Figure 3, the course of the NN content was significantly affected by the interaction of NFS and N rates. For the M-NFS, the data obtained fitted the quadratic regression model the best, characterized by an N optimum of 176 kg ha −1 and the respective maximum N nitrate content at BBCH 30 (N30max) of 46.7 kg ha −1 . This model clearly demonstrates that at the onset of FL, resources of NN were exploited more strongly, exceeding the N rate of 176 kg ha −1 . A quite different course of NN was observed for the O-NFS and OM-NFS treatments. The NN content increased progressively with the increased N rate. As compared to the quadratic model, characterizing the M-NFS, the linear model as obtained for the O-NFS shows nitrate N resources to be still present in the soil at the onset of FL in plots fertilized with highest N rates.
The difference in the course of the NN content during STME is best described by its quantitative change. As a rule, in both years, the NN content during this period significantly decreased. The effect of N rates was significant in both years, and the NN content decreased progressively with the increasing N rates. The NN content at BBCH 60, i.e., at the onset of WOSR flowering (FL), was significantly driven by both experimental factors and the interaction between them. It is necessary to stress that the NN content at end of the STME was significantly correlated with its status at the beginning of this phase, i.e., BBCH 30 (Table A1). In 2016, the impact of NFS on the NN content, evaluated independently of N rates, was the same as at BBCH 30. In 2017, a significantly higher amount of NN was recorded in the M-NFS. The effect of the progressively increased N rates was significantly different in both years. In 2016, the NN content increased up to the N rate of 120 kg ha −1 and then fell significantly, whereas in 2017 it increased in accordance with the increased N rates. As shown in Figure 3, the course of the NN content was significantly affected by the interaction of NFS and N rates. For the M-NFS, the data obtained fitted the quadratic regression model the best, characterized by an N optimum of 176 kg ha −1 and the respective maximum N nitrate content at BBCH 30 (N 30max ) of 46.7 kg ha −1 . This model clearly demonstrates that at the onset of FL, resources of NN were exploited more strongly, exceeding the N rate of 176 kg ha −1 . A quite different course of NN was observed for the O-NFS and OM-NFS treatments. The NN content increased progressively with the increased N rate. As compared to the quadratic model, characterizing the M-NFS, the linear model as obtained for the O-NFS shows nitrate N resources to be still present in the soil at the onset of FL in plots fertilized with highest N rates.  The dominant trends were evaluated for the NFS x N interaction ( Figure 4). For the M-NFS, the cubic regression model fitted the obtained data the best. An initial drop in the NN content of 55 kg ha −1 on the AC plot underwent stabilization at this level on plots fertilized with an N rate extending from 60 to 180 kg ha −1 . A marked decrease in the NN content was only recorded on the plot with 240 kg ha −1 of applied N. The course of the cubic model obtained for the O-NFS was completely different and can be divided into three stages. The first phase covers the stabilized range of the NN content. It amounted to −20 kg ha −1 of NN. In the next phase, NN content showed a progressive drop, which lasted up to the N rate of 175 kg ha −1 , reaching the lowest value of −106 kg ha −1 . In the third stage, resulting from the application of an N rate of 240 kg ha −1 , the NN content increased again, and the final decrease in its content with respect to onset of the STME was only −65 kg ha −1 . The trend obtained fully explains the progressive response of the NN content at the onset of FL to the progressively increased N fertilizer rates. The trend in the NN content for the OM-FS followed the quadratic regression model. The stability phase at the level of 16 kg ha −1 lasted only to the theoretical N rate of 29.3 kg ha −1 . Following this N value, an exponential decrease in the NN content along with increasing N rates was recorded. The difference in the course of the NN content during STME is best described by its quantitative change. As a rule, in both years, the NN content during this period significantly decreased. The effect of N rates was significant in both years, and the NN content decreased progressively with the increasing N rates.
The dominant trends were evaluated for the NFS x N interaction ( Figure 4). For the M-NFS, the cubic regression model fitted the obtained data the best. An initial drop in the NN content of 55 kg ha −1 on the AC plot underwent stabilization at this level on plots fertilized with an N rate extending from 60 to 180 kg ha −1 . A marked decrease in the NN content was only recorded on the plot with 240 kg ha −1 of applied N. The course of the cubic model obtained for the O-NFS was completely different and can be divided into three stages. The first phase covers the stabilized range of the NN content. It amounted to −20 kg ha −1 of NN. In the next phase, NN content showed a progressive drop, which lasted up to the N rate of 175 kg ha −1 , reaching the lowest value of −106 kg ha −1 . In the third stage, resulting from the application of an N rate of 240 kg ha −1 , the NN content increased again, and the final decrease in its content with respect to onset of the STME was only −65 kg ha −1 . The trend obtained fully explains the progressive response of the NN content at the onset of FL to the progressively increased N fertilizer rates. The trend in the NN content for the OM-FS followed the quadratic regression model. The stability phase at the level of 16 kg ha −1 lasted only to the theoretical N rate of 29.3 kg ha −1 . Following this N value, an exponential decrease in the NN content along with increasing N rates was recorded.
The third studied cardinal phase of WOSR growth with respect to the NN content was ripening (RE). The NN status at the RE phase was significantly correlated with its status at both the rosette and the onset of flowering stages, respectively (Table A1). It can be therefore concluded that the general pattern of NN management by WOSR did not change significantly from the onset of the STME phase. The course of NN content with respect to the increasing N rates was described by different regression models ( Figure 5). For the M-NFS, the linear regression model showing a progressive increase in the NN content fitted the data obtained the best. The cubic regression model best described treatments fertilized with organic N, but its course was treatment specific. For the O-NFS, three distinct phases of NN content were distinguished. The first phase, stabilized at the level of 16 kg ha −1 of NN, was related to an N rate of 15 kg ha −1 . The second phase, with a progressive increase in NN content, extended up to an N rate of 190 kg ha −1 . The third phase, characterized by an NN content decrease, covered the part of the curve with the highest N fertilizer rates. An opposite trend was found for the OM-NFS. The NN content, in spite of a slight variability, showed up to the N rate of 150 kg ha −1 a stabilization at the level of 30 kg ha −1 . A sudden increase in NN content was observed on the plot with the highest N rate.  The third studied cardinal phase of WOSR growth with respect to the NN content was ripening (RE). The NN status at the RE phase was significantly correlated with its status at both the rosette and the onset of flowering stages, respectively (Table A1). It can be therefore concluded that the general pattern of NN management by WOSR did not change significantly from the onset of the STME phase. The course of NN content with respect to the increasing N rates was described by different regression models ( Figure 5). For the M-NFS, the linear regression model showing a progressive increase in the NN content fitted the data obtained the best. The cubic regression model best described treatments fertilized with organic N, but its course was treatment specific. For the O-NFS, three distinct phases of NN content were distinguished. The first phase, stabilized at the level of 16 kg ha −1 of NN, was related to an N rate of 15 kg ha −1 . The second phase, with a progressive increase in NN content, extended up to an N rate of 190 kg ha −1 . The third phase, characterized by an NN content decrease, covered the part of the curve with the highest N fertilizer rates. An opposite trend was found for the OM-NFS. The NN content, in spite of a slight variability, showed up to the N rate of 150 kg ha −1 a stabilization at the level of 30 kg ha −1 . A sudden increase in NN content was observed on the plot with the highest N rate. The observed differences, both in seed yield and NN content at BBCH 30 and during the STME in response to the tested treatments, can be explained by analysis of the partial factor productivity of N. In this study, it was developed based on the amount of NN in a soil layer of 0.9 m at the onset of STME. In both years, the index values averaged over N rates were slightly higher for O-NFS and OM- The observed differences, both in seed yield and NN content at BBCH 30 and during the STME in response to the tested treatments, can be explained by analysis of the partial factor productivity of N. In this study, it was developed based on the amount of NN in a soil layer of 0.9 m at the onset of STME. In both years, the index values averaged over N rates were slightly higher for O-NFS and OM-NFS treatments. This dependence was significant for the OM-NFS treatment in 2017. The impact of N rates, averaged over NFSs, was significant for both years. The impact of NFS on the partial factor productivity of N-NO 3 (PFP N30 ) indices is shown in Figure 6. The key differences between NFSs result from the index value on the plot fertilized with an N rate of 60 kg ha −1 . In the M-NFS, the index was stable on this plot with respect to the N control. On treatments fertilized with organic N, the index increased by 31% for the O-NFS, and by 71% for the OM-FS. Following an N rate of 60 kg ha −1 , PFP N30 indices fell in response to increasing N rates, showing strong differences between NFSs, as proved by the regression models developed: Agronomy 2020, 10, x 11 of 29 Figure 6. Partial factor productivity of nitrate N at BBCH 30 on the background of nitrogen fertilization systems (NFSs). M-NFS, O-NFS, OM-NFS-mineral, organic, organic-mineral, respectively. a, b, c, d, e, f, g significance letters, a-the highest, g-the lowest; a means within a column followed by the same letter indicate a lack of significant difference between the treatments.

The Onset of Stem Elongation (STME)
The general rule governing the content of available P and K was its decline with soil depth. The content of available P in the topsoil at the onset of STME was, in general, in the suitable class (Table  4, [41]). In both years, a significant impact of NFSs on the P content was recorded in the first subsoil layer (0.3-0.6 m). Significantly lower values were recorded in treatments fertilized with organic N (O-NFS, OM-NFS). The effect of N rates was significant for the third soil layer in 2016, and for both subsoil layers in 2017. The lowest content of P was, as a rule, recorded in treatments fertilized with moderate N rates (120 and 180 kg N ha −1 in 2017 and 120 kg N ha −1 in 2017), and the highest in the N control plot. The content of available K in the topsoil was in the suitable class in 2016 and in the low in 2017. In the first year of the study, a significant impact of NFSs was observed for the third soil layer, and in the second in the whole subsoil. The general rule was the same as that observed for P, i.e., a significantly higher content of K was recorded in the M-NFS, and the effect of N rates was insignificant. Figure 6. Partial factor productivity of nitrate N at BBCH 30 on the background of nitrogen fertilization systems (NFSs). M-NFS, O-NFS, OM-NFS-mineral, organic, organic-mineral, respectively. a, b, c, d, e, f, g significance letters, a-the highest, g-the lowest; a means within a column followed by the same letter indicate a lack of significant difference between the treatments.
The PFP N30 models for treatments fertilized with organic N, in contrast to M-NFS, showed a stabilization phase, as indicated by the fixed PFP N30 minimum. For the O-NFS, its lowest value of 30.6 kg seeds kg −1 N was achieved for an N rate of 185.3 kg ha −1 . For the OM-NFS, these characteristics were 25.6 kg seeds per kg −1 N and 223.7 kg ha −1 , respectively. The PFP N30 index was, as a rule, negatively correlated with the NN content in all stages of WOSR growth (Table A1).

The Onset of Stem Elongation (STME)
The general rule governing the content of available P and K was its decline with soil depth. The content of available P in the topsoil at the onset of STME was, in general, in the suitable class (Table 4, [41]). In both years, a significant impact of NFSs on the P content was recorded in the first subsoil layer (0.3-0.6 m). Significantly lower values were recorded in treatments fertilized with organic N (O-NFS, OM-NFS). The effect of N rates was significant for the third soil layer in 2016, and for both subsoil layers in 2017. The lowest content of P was, as a rule, recorded in treatments fertilized with moderate N rates (120 and 180 kg N ha −1 in 2017 and 120 kg N ha −1 in 2017), and the highest in the N control plot. The content of available K in the topsoil was in the suitable class in 2016 and in the low in 2017. In the first year of the study, a significant impact of NFSs was observed for the third soil layer, and in the second in the whole subsoil. The general rule was the same as that observed for P, i.e., a significantly higher content of K was recorded in the M-NFS, and the effect of N rates was insignificant.
The content of available Mg in the topsoil was in the low class in 2016 and did not show any response to the experimental treatments. It is necessary to stress that its content, in contrast to P, increased with soil depth. In 2017, a significant impact of NFSs was recorded, but only for the third soil layer. A significantly higher Mg content was recorded in treatments fertilized with organic N, moving these two treatments to the suitable class of Mg availability. The effect of N rates was insignificant. The content of available calcium (Ca) was, in general, low (extremely low) and responded to applied treatments, but only in the wet 2017. The highest impact of NFSs was recorded in the topsoil. A significant increase with respect to the M-NFS was recorded on plots fertilized with organic N. The effect of N rates was observed in each soil layer. As a rule, the highest Ca content was recorded in soil fertilized with an N rate of 120 and 180 kg ha −1 .
For the whole N fertilization system, analyzed irrespectively of the N source, the first three Principal Components (PCs) explained 68.6% of the total variance. PC1 was significantly associated with four of fifteen variables. The highest positive loadings were recorded for Ca and Mg in the subsoil (Ca b , Ca c , Mg b ). A significant but negative loading was found for the K content in the second subsoil layer (K c ) (Table A2). PC2 had high, but negative loadings with P and Mg contents in the topsoil (P a , Mg a ) and with P in the first subsoil layer (P b ). PC3 was significantly but negatively associated with the nitrate N content at BBCH 30 N 30 and yield (Y). A significant impact on yield was exerted by K c and N 30 (Table S1a; Figure S1a). The applied stepwise analysis showed however a significant dependence of Y on Mg c , and N 30 : The four PCs explained 90.6% of the total variance in the M-NFS. PC1 was significantly correlated with six of fifteen analyzed variables (Table A2). Positive coefficients of correlation (for coefficient of determination -R 2 > 0.50) were recorded for Ca a,b,c , yield, and negative for P c , K c . PC2 was significantly but negatively correlated with P a , P b , Mg a . PC3 was significantly and positively associated with N 30 , and PC4 with Mg c . The strongest negative correlation with PC1 and Y was exerted by K in the third soil layer (K c ). The third group of variables associated with PC2 were contents of available P, Mg a , and N 30 . The fourth group of variables was represented by Mg c and PFP N30 showed another direction with respect to the previous group (Table S1b; Figure S1b).
The three PCs accounted for 82.2% of the total variance for the O-NFS. PC1 had significant loadings with six of fifteen variables. The contents of Ca b and Mg (Mg a , Mg b ) followed the same direction. PC2 had positive loadings with Mg c and PFP N30 , and negative with N 30 . PC3 was associated with Y, showing a significant relationship with N 30 , but at the same time negative with P c (Table S1c; Figure S1c). The three PCs explained 85.9% of the OM-NFS total variance. PC1 was significantly and positively associated with Mg b , Mg c , and Ca for all soil layers (Table S1d; Figure S1d). Negative loadings were exerted by K b , K c . PC2 had high loadings with P a , P b and Mg a . PC3 was associated with N 30 and Y, but they had negative loadings. The highest positive impact on Y was exerted by N 30 .

The Onset of Flowering-Change in Nutrient Availability during STME
The analysis of soil fertility at the onset of FL comprises two components: (i) current status of the content of available nutrient, (ii) the net content of available nutrient change during STME (Tables 4 and 5). In both years, the content of available P in the topsoil was not affected by the studied NFSs. As compared to the onset of STME, the P content increased. The effect of NFSs was observed in both years, and a significant increase due to application of organic N was recorded in 2017. The effect of N rates was significant in 2017. The P content in the subsoil underwent a great change during STME. In 2016, it was depleted, being the strongest in the M-NFS. In 2017, an opposite trend was recorded. As a result, in 2016, a significantly higher P content at the onset of FL was recorded in treatments with organic N, whereas in 2017, with mineral N. In 2016, the lowest P depletion with respect to the N rates was recorded on the plot fertilized with 120 kg ha −1 of N. In 2017, the net P content increase was variable, but the highest was recorded in the plot with 240 kg ha −1 of N, leading to its highest content at the onset of FL.
The content of available K showed a much lower variability during STME as compared to P. In both years, its highest content was recorded in the M-NFS. The significant differences between NFSs as observed in 2016 were due to an extremely high increase in the K content in soil fertilized with mineral N. In 2017, a net K content increase with soil depth was recorded only in the O-NFS treatment. The application of mineral N resulted in K depletion in the topsoil, which was compensated by its increase in the subsoil, especially in the OM-NFS. The effect of N rates on both the K content and its change during STME was not significant. However, in 2017, a year with a significantly lower initial K content, a decreasing trend in accordance with the increasing N rates was observed in the topsoil. The K status at the onset of FL was a result of the decreasing trend in the net K content increase during STME, which underwent depletion in plots fertilized with the highest N rate.
At the onset of FL, the content of available Mg showed a great change with respect to the onset of STME and responded significantly to NFSs in the subsoil. In 2016, as in the case of K, a significantly higher Mg content was recorded in the M-NFS. The key reason for the observed Mg status was a net increase in its content, especially in the deepest soil layer. In 2017, an opposite trend was observed, i.e., the content of available Mg was depleted, being the strongest in the M-NFS. The increasing N rates significantly affected the content of available Mg in the third soil layer in 2016, and in both subsoil layers in 2017. In both years, the content of available Ca showed an extremely high variability in response to the experimental factors. In 2016, in the topsoil and in the first subsoil layer, a significantly higher content of Ca was recorded in the M-NFS, whereas in 2017, in the OM-NFS. The effect of N rates on the content of Ca in both years was recorded for the topsoil. In 2016, the observed difference between N plots was due to significant differences in the net Ca content increase, as occurred in plots fertilized with 60, 120, and especially with 240 kg ha −1 of N. In 2017, the Ca content at the onset of FL resulted from a great variability in its content change during STME. A net increase was the attribute of plots fertilized with an N rate of 60 and 240 kg ha −1 . In the other plots, the Ca content underwent depletion. In both years, these trends were also observed in the deeper soil layers.
The pooled analysis of NFSs show that three PCs accounted for 58.64% of the total variance (Table A3). PC1 was associated with five of fifteen variables, of which Mg had positive loadings, irrespective of soil layer. P showed the opposite direction on the PC1 axis in the subsoil. PC2 was significantly associated with Ca c , and PC3 with Y. The impact of soil variables on Y was insignificant, and the highest positive was exerted only by the nitrate N content at BBCH 60 (N 60 ), which in turn was positively affected by Mg (Table S2a; Figure S2a). For the M-NFS, two PCs were extracted that cumulatively explained 59.35% of the total variance. PC1 accounted for 45.2% of the total variance and had high positive loadings for Mg, including all soil layers, and for K, but only for the topsoil (K a ). Negative loadings were recorded for P, also including all soil layers. PC2 was not associated with any variable with high loading (R 2 > 0.50). Yield was weakly related to the soil variables but N 60 showed a significant and positive response to Mg in the subsoil and Ca in the second subsoil layer (Table S2b; Figure S2b).
Three PCs explained 73.24% of the total variance for the O-NFS (Table A3). PC1 was dominated by six of fifteen variables that had high loadings. A positive impact on PC1 was exerted only by Y, and a negative one by PFP N30 , and also by Mg and Ca in the subsoil. The second set of variables changed in the opposite direction to Y. PC2 was associated with P (P a , P b ) that had high and positive loadings. PC3 grouped two variables, i.e., K a with a negative and Ca a with a positive loading. The efficiency of N 30 was significantly affected by Ca and Mg contents in the second subsoil layer. Yield was positively, but not significantly, dependent on N 60 , but negatively on the contents of Ca and Mg in the first subsoil layer (Table S2c; Figure S2c).
Four PCs accounted for 85.07% of the total variance for the OM-NFS. PC1 was associated with six of fifteen variables, such as K a , K c , and Mg a , which had high and positive loadings. A negative loading was exerted by Ca b . PC2 showed high and positive loadings with Mg b and P b . PC3 had negative loadings with Ca a and Y. PC4 was significantly, but negatively, associated with P a . The NN content at the onset of FL was significantly correlated with Mg a , and negatively with P b . The content of Ca in the topsoil and N 60 showed the highest and at the same time positive relationships with yield (Table S2d; Figure S2d).

Physiological Maturity
The content of available P at the physiological maturity of WOSR as compared to the onset of FL, as a rule, was slightly higher in 2016, and considerably lower in 2017. A significant effect of NFSs was observed only in 2016 (Table 6). A higher P content was recorded for treatments fertilized with organic N, being significant in the subsoil. The effect of N rates was non-significant, although a decreasing trend in accordance with the increasing N rates was observed. The content of available K was only affected by NFS. In 2016, as compared to FL, much lower K values were recorded in the topsoil, but an opposite trend was observed in the subsoil. In 2017, an increase in K content was noted for the subsoil.
A depletion of the available Mg content, as compared to FL, was recorded in 2016 for all soil layers. In 2017, it appeared only in the topsoil. In the subsoil, a net increase in Mg content was recorded. The impact of NFS was observed only in 2016 for the subsoil, in which a significantly higher content of Mg was recorded for the M-NFS and O-NFS. The effect of N rates was observed only for the second subsoil layer in 2016. The content of available Ca showed the highest variability. Ca status was significantly affected by both experimental factors and their interaction with a particular year. In both years, the impact of NFS was recorded for all three soil layers. As a rule, the lowest content of available Ca was found for the OM-NFS. The highest Ca values were the attribute of M-NFS. The effect of the increasing N rates was highly variable between soil layers. For the topsoil, the lowest content of Ca was generally recorded for the plot fertilized with 180 kg ha −1 of N. In the subsoil, no consistent rule governing the content of available Ca was observed.
The four PCs accounted for 69.50% of the total variance in the pooled NFSs (Table A4). PC1 had the highest and most positive loadings for the Ca content in the subsoil. PC2 was associated with K c that exerted a negative effect in its value. PC3 explained 14.16% of the total variance, and no significant loading was recorded. PC4 was controlled by the nitrate N content at BBCH 89 (N 89 ), having a negative loading. Yield showed a positive relationship with N 89 and P a , but the latter changed in the opposite direction to that of N 89 (Table S3a; Figure S3a). For the M-NFS, three PCs accounted for 75% of the total variance. PC1 was associated with seven of fourteen variables. Positive loadings were identified for P in the subsoil and Ca a , but negative values for Mg and Ca for the subsoil. Variables within each subgroup were significantly and positively correlated with each other. Yield exerted a high positive loading on PC2. PC3 had the highest and most positive loading with K a , and K b . Yield was weakly related to other soil variables but showed an opposite direction to P a and K c . The latter variable was significantly and negatively correlated with Y (Table S3b; Figure S3b).
The O-NFS was associated with three PCs, which explained 73.5% of the total variance. PC1 showed high loadings of P b , which were positive, and of Mg b , Ca b , Ca c , which were negative. PC2 showed high positive loading with N 89 and moderate with Y (0.65) and Kc (0.63). These three variables were significantly correlated with each other, although only N 89 significantly affected Y. The opposite direction to that set of variables was exerted by P variables, of which P a and P b were significantly and negatively correlated with Y. PC3 showed a high but negative loading with K a (Table S3c; Figure S3c).
In the OM-NFS, four PCs accounted for 81.88% of the total variance. PC1 had positive loadings for Ca, but negative for P and Mg b in the subsoil. PC2 was significantly affected by Mg a and Mg c , which had high and positive loadings, and P a with a negative effect. PC3 was dominated with K a . PC4 was associated with N 89 due to its high loading and moderately affected by Y (0.61). Yield was closely, but not significantly, related to K a and N 89 (Table S3d; Figure S3d).

In-Season Variability in the Growth Factor-Nitrate Nitrogen
The yield variability of winter oilseed rape (WOSR) was evaluated based on the in-season variability in the N-NO 3 content, i.e., the nitrogenous growth factor (NGF), directly impacting the rate of WOSR plant growth during the growing season. Yield responded to the applied N, irrespective of its chemical form, but significantly only to its lowest rate of 60 kg ha −1 . This N rate resulted in a yield increase of 60% in 2016, a year with drought in June, and a 115% increase in 2017. The recorded yield increase fully corroborates the opinion on the high sensitivity of WOSR to drought during the end of inflorescence emergence and the onset of flowering [18,39]. The strong response of WOSR to the lowest N rate is typical for soil naturally poor in inorganic N, or soils strongly depleted of this N form [43]). The second problem, which became apparent during the study, refers to the lack of WOSR response to increasing N rates, applied as ammonium nitrate (AN) or digestate, which in the basic form is ammonium [44]. The application of 240 kg ha −1 of N as AN resulted in a yield reduction. The observed tendencies are broadly explained by assuming an imbalance of applied nutrients or referring to the low nitrogen use efficiency (NUE) of WOSR plants [10,[45][46][47][48][49]. In the studied case, this observation was fully corroborated by the trend of partial factor productivity of N indices, which were calculated based on the N-NO 3 content in the soil at BBBCH 30. The observed PFP N30 stabilization on treatments with the highest N rates was due to the lower rate of the NN content decrease during STME, subsequently leading to the better N utilization by WOSR during the seed filling period (SFP). The net N increase during the SFP, as has been documented recently by Łukowiak and Grzebisz [43], is the prerequisite of the high yield of WOSR.
One of the most important targets of soil or plant N monitoring during the growing season is to make a reliable yield prognosis [3,50]. In the study, the N-NO 3 content in the soil profile during two of the cardinal stages of WOSR growth, i.e., rosette (BBCH 30) and the onset of flowering (BBCH 60), significantly limited seed yield, and therefore can be used as a yield predictor (Table A1). The relationship obtained was best described using the quadratic regression model: As results from equations developed, the N-NO 3 content of 150 kg ha −1 , as recorded at BBCH 30, resulted in a theoretical Y max of 3.91 t ha −1 . The calculated yield is very close to that achieved in the OM-NFS on the plot with the N rate of 180 kg ha −1 . The yield prognosis at the onset of flowering was significantly lower when compared to BBCH 30, as indicated by the Y max of 3.451 t ha −1 . This yield could be achieved provided the N-NO 3 amounted to 47.7 kg ha −1 . Based on these two data sets, it can be concluded that a typical feature of N management by WOSR during the stem elongation phase is a significant decline in the N-NO 3 content. The average N-NO 3 content decrease, based on the respective N optima, was 102.3 kg ha −1 (= 150 kg ha −1 at BBCH 30 minus 47.7 kg ha −1 at BBCH 60). The study showed that the change in the N-NO 3 content (∆N-NO3) within the STME can be determined based on the N-NO 3 content precisely at BBCH 30: The direction coefficient of the linear regression model obtained clearly indicates a high efficiency of N-NO 3 during STME, which reached 85% of the N-NO 3 content in the whole analyzed soil profile. This value can be used as an index of nitrate N uptake efficiency by WOSR.
The analysis of the N-NO 3 content at BBCH 30 and its rate of change during STME can be used as a basis for the revision of the current hypothesis about the critical stages of WOSR yield development. In fact, it is well-documented that two phases, i.e., inflorescence emergence and flowering, are crucial for the establishment of basic yield components. Any disturbance in plant growth during these two phases, either abiotic or nutritional, leads to yield reduction [18,39,51]. The observed dependence of the N-NO 3 decrease during STME on its content in the soil profile at BBCH 30, as documented in this study, explicitly indicates the rosette stage as the cardinal for exploiting WOSR yielding potential [52]. In the well-established WOSR canopy, as in this study, the NN content underwent significant depletion, irrespective on its amount and N form (Figure 7). OM-NFS on the plot with the N rate of 180 kg ha −1 . The yield prognosis at the onset of flowering was significantly lower when compared to BBCH 30, as indicated by the Ymax of 3.451 t ha −1 . This yield could be achieved provided the N-NO3 amounted to 47.7 kg ha −1 . Based on these two data sets, it can be concluded that a typical feature of N management by WOSR during the stem elongation phase is a significant decline in the N-NO3 content. The average N-NO3 content decrease, based on the respective N optima, was 102.3 kg ha −1 (= 150 kg ha −1 at BBCH 30 minus 47.7 kg ha −1 at BBCH 60). The study showed that the change in the N-NO3 content (ΔN-NO3) within the STME can be determined based on the N-NO3 content precisely at BBCH 30: ΔN-NO3 = −0.85N30 + 23.4 for n = 15, R 2 = 0.98, P ≤ 0.01 (11) The direction coefficient of the linear regression model obtained clearly indicates a high efficiency of N-NO3 during STME, which reached 85% of the N-NO3 content in the whole analyzed soil profile. This value can be used as an index of nitrate N uptake efficiency by WOSR.
The analysis of the N-NO3 content at BBCH 30 and its rate of change during STME can be used as a basis for the revision of the current hypothesis about the critical stages of WOSR yield development. In fact, it is well-documented that two phases, i.e., inflorescence emergence and flowering, are crucial for the establishment of basic yield components. Any disturbance in plant growth during these two phases, either abiotic or nutritional, leads to yield reduction [18,39,51]. The observed dependence of the N-NO3 decrease during STME on its content in the soil profile at BBCH 30, as documented in this study, explicitly indicates the rosette stage as the cardinal for exploiting WOSR yielding potential [52]. In the well-established WOSR canopy, as in this study, the NN content underwent significant depletion, irrespective on its amount and N form (Figure 7).

In-Season Variability in Soil Fertility Factor-Available Nutrients
The yield of WOSR as discussed in the previous section responded significantly to the applied fertilizer N, irrespective of its chemical form. The key question remaining is to what the extent the soil fertility factor (FF), as defined by the in-season status of four basic nutrients (P, K, Mg, Ca), was related to the content of N-NO3, defining the NGF status, and in consequence, seed yield. The applied PCA was revealed as a useful tool to evaluate the nitrogen fertilization system based on the above described factors. It was also a useful tool to define the production role of particular subsoil layers with respect to the availability of a given nutrient, or set of nutrients, as yield limiting factors. So far,

In-Season Variability in Soil Fertility Factor-Available Nutrients
The yield of WOSR as discussed in the previous section responded significantly to the applied fertilizer N, irrespective of its chemical form. The key question remaining is to what the extent the soil fertility factor (FF), as defined by the in-season status of four basic nutrients (P, K, Mg, Ca), was related to the content of N-NO 3 , defining the NGF status, and in consequence, seed yield. The applied PCA was revealed as a useful tool to evaluate the nitrogen fertilization system based on the above described factors. It was also a useful tool to define the production role of particular subsoil layers with respect to the availability of a given nutrient, or set of nutrients, as yield limiting factors. So far, knowledge about nutrient availability from deeper soil layers, especially during the growing season, is low [53].
The analysis of the NGF, i.e., variability in the N-NO 3 content during the growing season, clearly indicates the rosette and the onset of flowering stages as cardinal for determining WOSR yield. Based on the PCA, it can be concluded that at the onset of STME, the primary set of soil nutrients, including the content of available Ca and Mg, representing the FF significantly limited yield, irrespective of the studied NFS. It is necessary to stress that the importance of both nutrients for NFS productivity resulted more strongly from their contents in the subsoil than in the topsoil. This conclusion stresses the importance of Ca as the yield nutritional factor. In fact, the yield forming function of Ca, in spite of the high requirements of some crops, such as oilseed rape, is weakly recognized [23,53]. In contrast to Ca, the yield forming functions of Mg are well recognized, but are mostly related to the onset of flowering and during the SFP [18,49].
The detailed analysis of the impact of FFs on yield for each of the studied NFSs showed a significantly different response to the applied N form. The general pattern of the FF relation with Y, as presented above, was to a great extent typical for the M-NFS. The shortage of Ca and Mg in connection with a shortage of NN resulted in the insufficient exploitation of both P and especially K from the second subsoil layer. The content of available K in this layer in the M-NFS can be used as a single yield (Y) predictor at BBCH 30: Y = −0.68K c + 6.32 for n = 10, R 2 = 0.85, P ≤ 0.01 (12) The negative sign of the direction coefficient implicitly indicates that the low exploitation of K from the subsoil, just at the rosette stage of WOSR, is the reason for the yield decline. The OM-NFS yielded higher, because N was better balanced with basic nutrients, such as P and K, which were more strongly exploited by WOSR plants from the subsoil as compared to the M-NFS (Table S1b and Table 3d). It was found that a higher content of available K in the subsoil resulted in favorable conditions for the uptake of Ca and Mg. This hypothesis was fully corroborated at the onset of flowering, when available Ca positively affected N availability, finally resulting in a yield increase. The same effect in the OM-NFS on N availability was exerted by P and Mg. It is well-documented that all these three nutrients are strongly exploited by WOSR from subsoil [54][55][56][57].
A significantly different pattern of FF impact on the content of NN and yield was observed in the O-NFS. In this particular NFS, yield significantly depended on the content of N-NO 3 , i.e., the NGF was found to be the direct yield driver (Table S1c; Figure S1c). The productivity of this system was also dependent on the content of Mg and Ca in the deep subsoil layer (c). The higher availability of both nutrients in this layer significantly affected the unit productivity of N (PFP N30 ). This was the key reason for the much higher productivity of the O-NFS system as compared to the M-NFS.
This study explicitly shows that a reasonable, but significant, decrease in the NN content during STME is the prerequisite of high WOSR yield (Figure 7). A similar decrease in the content of other nutrients can therefore be assumed. K requires special attention, because during this period it reaches both the maximum rate of uptake and the maximum value of its accumulation just at the end of the inflorescence emergence (INFE) [23,56]. In 2016, a net K content increase was recorded in all three soil layers. The same trend was observed in 2017 for treatments fertilized fully or partly with organic N. The observed phenomenon was not, however, related to the amount of K applied in digestate. Therefore, it cannot be explained by the direct impact of this fertilizer on K availability through K input, or indirectly by ammonium oxidation and its subsequent impact on the cation exchange [58]. The only reasonable explanation for the K content increase with soil depth is the WOSR activity in the rhizosphere, including its acidification [59].
At the onset of flowering, being at the same time the end phase of intensive uptake of N by WOSR, the key limiting nutrient for the whole system became the content of available Mg, which controlled the amount of NN, the direct yield driver. The result obtained corroborates the importance of Mg for stabilizing the nutritional status of WOSR plants at this particular phase [18]. The insufficient availability of Mg on the one hand and excess of available P in the subsoil on the other points to FFs as critical factors for WOSR yield performance, irrespective of the NFS. Magnesium content and to a lesser extent Ca in the subsoil exerted a positive impact on the content of NN, a direct growth factor. Therefore, it can be concluded that the higher the content of available Mg in the subsoil, the better the supply of N to WOSR plants during the INFE, the phase responsible for seed set [49]. In the O-NFS, the impact of the key FF variables on the NFS functioning revealed a positive impact of K c on the content of NN at the end of the INFE phase, concomitant with a significant improvement in N productivity, which was exerted by the Mg and Ca present in the second subsoil layer ( Table S2c). The positive impact of Mg content on the NN content was fully corroborated in the OM-NFS. The higher content of both nutrients at the onset of flowering resulted in a positive yield response. The yield increase also resulted from the concomitant positive impact of available Ca on N productivity. In both years, the highest net increase in the content of available Ca during the STME was recorded on plots fertilized with the highest N rate, which resulted in the highest yield. These three simultaneously occurring processes led to a higher seed yield and the advantage of the OM-NFS over other tested NFSs.
The dominant impact of FFs as limiting production factors was fully corroborated at WOSR physiological maturity. The entire NFS production efficiency was governed by the content of available Ca and P in the subsoil. It is well documented that in soil of a neutral pH range, the content of available P depends on the content of Ca and vice versa (Table S3a) [17]. The negative relationship between the content of available K in the second subsoil layer with yield clearly indicates the stability of the soil/plant system during the spring growing season: Y = −0.0015K c 2 + 0.15K-0.334 for n = 10, R 2 = 0.55, P ≤ 0.05 (13) The presented equation shows that the depletion of the K content in the second subsoil layer at WOSR physiological maturity to 50.3 mg kg −1 is the prerequisite of the maximum yield of 3.466 t ha −1 . A quite different set of variables were responsible for O-NFS. In the final phase of WOSR growth, the two decisive variables for yield were the amount of NN, K, and also P. The importance of this set of nutrients for exploiting WOSR yielding potential corroborates the latest study by Grzebisz et al. [49].
The key importance of K during the SFP results from its physiological function, as related to the transport of assimilates to the growing pods and seeds [60]. The better the supply of assimilates to the growing pods and seeds, the higher the number of seeds, subsequently resulting in a higher yield. The shortage of NN at this stage led to the weak exploitation of P from the whole soil profile. This observation is in agreement with Grzebisz et al. [61], who showed that the exploitation of P from its soil resources depends on the WOSR sink strength, which is related to the number of seeds per unit area.

Conclusions
The rosette stage of WOSR growth (BBCH 30) was revealed as the cardinal for seed yield performance. The highest yield was obtained, provided two basic conditions were fulfilled. The first refers to the N-NO 3 content at this particular stage, amounting to 150 kg ha −1 . The second was the considerable decrease in the N-NO 3 content during the stem elongation phase of WOSR growth. The yield was significantly dependent on the nitrogen fertilization system and increased in the order of NFS: M, O, OM. The maximum yields of 3.616, 3.887, 4.195 t ha −1 were achieved given the optimum N rates of 138.9, 171.4, 210 kg ha −1 , respectively. The key factors limiting yield were significantly related to the particular NFS. In the case of M-NFS, yield was limited by the availability of Ca and Mg, but especially of K in the deepest subsoil layer (0.6-0.9 m). In the O-NFS, the key limiting factor was the shortage of N-NO 3 . The highest yield obtained in the OM-NFS was due to the positive impact of Ca and Mg on both the amount of N-NO 3 and its productivity. During the second cardinal stage of WOSR yield performance, i.e., at the onset of flowering, Mg and to a lesser extent the contents of Ca in the subsoil both resulted in an N-NO 3 content increase-a direct growth factor. In can be therefore concluded that the higher the content of available Mg in the subsoil at the onset of WOSR flowering, the better the supply and utilization of N by plants, resulting in the highest yield. PCA analysis was found to be an effective tool in discriminating against nitrogen fertilization systems differing in the source of N (mineral vs. organic).  Table S1a. Spearman's matrix of correlation coefficients for the pooled nitrogen fertilization systems at BBCH 30 of WOSR growth, n = 30; Table S1b. Spearman's matrix of correlation coefficients for the mineral nitrogen fertilization systems at BBCH 30 of WOSR growth, n = 10; Table S1c. Spearman's matrix of correlation coefficients for organic nitrogen fertilization systems at BBCH 30 of WOSR growth, n = 10; Table  S1d. Spearman's matrix of correlation coefficients for the organic-mineral nitrogen fertilization systems at BBCH 30 of WOSR growth, n = 10; Table S2a. Spearman's matrix of correlation coefficients for the pooled nitrogen fertilization systems at BBCH 60 of WOSR growth, n = 30; Table S2b. Spearman's matrix of correlation coefficients for the mineral nitrogen fertilization systems at BBCH 60 of WOSR growth, n = 10; Table S2c. Spearman's matrix of correlation coefficients for organic nitrogen fertilization systems at BBCH 60 of WOSR growth, n = 10; Table S2d. Spearman's matrix of correlation coefficients for the organic-mineral nitrogen fertilization systems at BBCH 60 of WOSR growth, n = 10; Table S3a. Spearman's matrix of correlation coefficients for the pooled nitrogen fertilization systems at BBCH 89 of WOSR growth, n = 30; Table S3b. Spearman's matrix of correlation coefficients for the mineral nitrogen fertilization systems at BBCH 89 of WOSR growth, n = 10; Table S3c. Spearman's matrix of correlation coefficients for organic nitrogen fertilization systems at BBCH 89 of WOSR growth, n = 10; Table S3d. Spearman's matrix of correlation coefficients for the organic-mineral nitrogen fertilization systems at BBCH 89 of WOSR growth, n = 10. Figure S1

Conflicts of Interest:
The authors declare no conflict of interest.