**3. Results**

#### *3.1. In-Season Variability in the Nitrate Nitrogen (N-NO3, 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 (<sup>1</sup>/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%.


**Table 3.** In-season variability of the nitrate nitrogen content variability.

\*\*\*, \*\*, \* significant at *p* < 0.001; < 0.01; < 0.05, respectively; n.s.—non significant; a, b, c, d significance letters, a—the highest, d—the lowest; a means within a column followed by the same letter indicate a lack of significant difference between the treatments. N30, N60, N89—nitrate N content at BBCH 30, 60, and 89, kg ha−<sup>1</sup> respectively; N30–N60—nitrate N balance at BBCH 60 vs. BBCH 30; PFPN30–partial factor productivity of nitrogen nitrate (N-NO3) at WOSR stage of BBCH 30 (rosette).

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−<sup>1</sup> did not result in any yield increase. An N rate of 240 kg ha−<sup>1</sup> 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:

$$\text{Y-M-NFS} = -0.00009 \text{N}^2 + 0.025 \text{N} + 1.88 \text{ for } n = 5, R^2 = 0.93, \text{P} \le 0.05 \tag{2}$$

$$\text{Y-O-NFS} = -0.00007 \text{N}^2 + 0.024 \text{N} + 1.83 \text{ for } n = 5, R^2 = 0.94, \text{P} \le 0.05 \tag{3}$$

$$\text{Y-OM-NFS} = -0.00005 \text{N}^2 + 0.021 \text{N} + 1.99 \text{ for } n = 5, R^2 = 0.97, \text{P} \le 0.05 \tag{4}$$

**Figure 1.** Yield of winter oilseed rape as affected by N rates on the background of nitrogen fertilization systems (NFSs). M-NFS, O-NFS, OM-NFS—mineral, organic, organic-mineral, respectively. a, b, c, d significance letters, a—the highest, d—the lowest; a means within a column followed by the same letter indicate a lack of significant difference between the treatments.

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−<sup>1</sup> was the attribute of the M-NFS. The respective Ymax amounted to 3.616 t ha−1. The parameters of the O-NFS were 171.4 kg ha−<sup>1</sup> of applied N, and the Ymax of 3.887 t ha−1. The highest values of both parameters were found for the OM-NFS. A Nop of 210 kg ha−<sup>1</sup> resulted in a theoretical Ymax 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−<sup>1</sup> 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−<sup>1</sup> 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.

**Figure 2.** Trends in the nitrate N content at the onset of WOSR stem elongation stage of growth for three nitrogen fertilization systems (NFSs). M-NFS, O-NFS, OM-NFS—mineral, organic, organic-mineral, respectively.

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−<sup>1</sup> 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−<sup>1</sup> 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.

**Figure 3.** Trends in the N nitrate content at the onset of WOSR flowering for three nitrogen fertilization systems (NFSs). M-NFS, O-NFS, OM-NFS—mineral, organic, organic-mineral NFS, respectively.

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−<sup>1</sup> 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−<sup>1</sup> 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−<sup>1</sup> 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−<sup>1</sup> 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 managemen<sup>t</sup> 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−<sup>1</sup> 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−<sup>1</sup> 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.

**Figure 4.** Trends in the N nitrate content change during the stem elongation phase of WOSR growth for three nitrogen fertilization systems (NFSs). M-NFS, O-NFS, OM-NFS—mineral, organic, organic-mineral, respectively.

**Figure 5.** Trends in the N nitrate content at WOSR ripening for three nitrogen fertilization systems (NFSs). M-NFS, O-NFS, OM-NFS—mineral, organic, organic-mineral NFS, respectively.

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-NO3 (PFPN30) 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, PFPN30 indices fell in response to increasing N rates, showing strong differences between NFSs, as proved by the regression models developed:

$$\text{M-NFS PFP}\_{\text{N50}} = -0.16 \text{N} + 57.7 \text{ for } n = 4, R^2 = 0.73, \text{P} \le 0.05 \tag{5}$$

$$\text{O-NFS: }\text{PFP}\_{\text{N30}} = 0.0017\text{N}^2 - 0.63\text{N} + 89 \text{ for } n = 4, R^2 = 0.99, \text{P} \le 0.01\tag{6}$$

$$\text{OM-NFS:}\;\text{PFP}\_{\text{NS}0} = 0.0019 \text{N}^2 - 0.85 \text{N} + 120.7 \text{ for } n = 4, \; \text{R}^2 = 0.92, \text{P} \le 0.01 \tag{7}$$

**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 PFPN30 models for treatments fertilized with organic N, in contrast to M-NFS, showed a stabilization phase, as indicated by the fixed PFPN30 minimum. For the O-NFS, its lowest value of 30.6 kg seeds kg−<sup>1</sup> 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−<sup>1</sup> N and 223.7 kg ha−1, respectively. The PFPN30 index was, as a rule, negatively correlated with the NN content in all stages of WOSR growth (Table A1).

#### *3.2. In-Season Variability in Contents of Available Nutrients*

#### 3.2.1. 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 e ffect 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−<sup>1</sup> in 2017 and 120 kg N ha−<sup>1</sup> 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 e ffect 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 e ffect 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 e ffect 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 (Cab, Cac, Mgb). A significant but negative loading was found for the K content in the second subsoil layer (Kc) (Table A2). PC2 had high, but negative loadings with P and Mg contents in the topsoil (Pa, Mga) and with P in the first subsoil layer (Pb). PC3 was significantly but negatively associated with the nitrate N content at BBCH 30 N30 and yield (Y). A significant impact on yield was exerted by Kc and N30 (Table S1a; Figure S1a). The applied stepwise analysis showed however a significant dependence of Y on Mgc, and N30:

$$\text{Y} = 1.852 - 0.025 \text{ Mg}\_{\text{c}} \text{ }^{\ast} + 0.012 \text{N}\_{\text{30}} \text{ }^{\ast \ast} \text{ for } n = 30 \text{, } R^2 = 86 \tag{8}$$

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 coe fficients of correlation (for coe fficient of determination - *R<sup>2</sup>* > 0.50) were recorded for Caa,b,c, yield, and negative for Pc, Kc. PC2 was significantly but negatively correlated with Pa, Pb, Mga. PC3 was significantly and positively associated with N30, and PC4 with Mgc. The strongest negative correlation with PC1 and Y was exerted by K in the third soil layer (Kc). The third group of variables associated with PC2 were contents of available P, Mga, and N30. The fourth group of variables was represented by Mgc and PFPN30 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 Cab and Mg (Mga, Mgb) followed the same direction. PC2 had positive loadings with Mgc and PFPN30, and negative with N30. PC3 was associated with Y, showing a significant relationship with N30, but at the same time negative with Pc (Table S1c; Figure S1c).


letter indicate a lack of significant difference between the treatments, a, b, c—soil layers of 0.0–0.30, 0.30–0.60, 0.60–0.90 cm, respectively.

**Table 4.** The content of available nutrients along the soil layers at the onset of stem elongation of winter oilseed rape (BBCH 30), mg kg−1 soil.

*Agronomy* **2020**, *10*, 1701

The three PCs explained 85.9% of the OM-NFS total variance. PC1 was significantly and positively associated with Mgb, Mgc, and Ca for all soil layers (Table S1d; Figure S1d). Negative loadings were exerted by Kb, Kc. PC2 had high loadings with Pa, Pb and Mga. PC3 was associated with N30 and Y, but they had negative loadings. The highest positive impact on Y was exerted by N30.

#### 3.2.2. 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 a ffected by the studied NFSs. As compared to the onset of STME, the P content increased. The e ffect of NFSs was observed in both years, and a significant increase due to application of organic N was recorded in 2017. The e ffect of N rates was significant in 2017. The P content in the subsoil underwent a grea<sup>t</sup> 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−<sup>1</sup> of N. In 2017, the net P content increase was variable, but the highest was recorded in the plot with 240 kg ha−<sup>1</sup> 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 di fferences 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 e ffect 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 grea<sup>t</sup> 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 a ffected 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 e ffect of N rates on the content of Ca in both years was recorded for the topsoil. In 2016, the observed di fference between N plots was due to significant di fferences in the net Ca content increase, as occurred in plots fertilized with 60, 120, and especially with 240 kg ha−<sup>1</sup> of N. In 2017, the Ca content at the onset of FL resulted from a grea<sup>t</sup> 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 Cac, 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 (N60), which in turn was positively a ffected by Mg (Table S2a; Figure S2a).


\*\*\*, \*\*, \* significant at *p* < 0.001; < 0.01; < 0.05, respectively; n.s.—non significant; a, b, c significance letters, a—the highest, c—the lowest; a means within a column followed by the same

letter indicate a lack of significant difference between the treatments, a, b, c—soil layers of 0.0–0.30, 0.30–0.60, 0.60–0.90 cm, respectively.

*Agronomy* **2020**, *10*, 1701

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 (Ka). Negative loadings were recorded for P, also including all soil layers. PC2 was not associated with any variable with high loading (*R<sup>2</sup>* > 0.50). Yield was weakly related to the soil variables but N60 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 PFPN30, 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 (Pa, Pb) that had high and positive loadings. PC3 grouped two variables, i.e., Ka with a negative and Caa with a positive loading. The efficiency of N30 was significantly affected by Ca and Mg contents in the second subsoil layer. Yield was positively, but not significantly, dependent on N60, 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 Ka, Kc, and Mga, which had high and positive loadings. A negative loading was exerted by Cab. PC2 showed high and positive loadings with Mgb and Pb. PC3 had negative loadings with Caa and Y. PC4 was significantly, but negatively, associated with Pa. The NN content at the onset of FL was significantly correlated with Mga, and negatively with Pb. The content of Ca in the topsoil and N60 showed the highest and at the same time positive relationships with yield (Table S2d; Figure S2d).
