Wheat Yield Estimation with NDVI Values Using a Proximal Sensing Tool

Nitrogen (N) splitting is critical to achieving high crop yields without having negative effects on the environment. Monitoring crop N status throughout the wheat growing season is key to finding the balance between crop N requirements and fertilizer needs. Three soft winter wheat fertilization trials under rainfed conditions in Mediterranean climate conditions were monitored with a RapidScan CS-45 (Holland Scientific, Lincoln, NE, USA) instrument to determine the normalized difference vegetation index (NDVI) values at the GS30, GS32, GS37, and GS65 growth stages. The threshold NDVI values in the Cezanne variety were 0.7–0.75 at the GS32, GS37, and GS65 growing stages. However, for the GS30 growing stage, a threshold value could not be established precisely. At this stage, N deficiency may not affect wheat yield, as long as the N status increases at GS32 stage and it is maintained thereafter. Following the NDVI dynamic throughout the growing season could help to predict the yields at harvest time. Therefore, the ΣNDVI from GS30 to GS65 explains about 80% of wheat yield variability. Therefore, a given yield could be achieved with different dynamics in wheat NDVI values throughout the growing cycle. The determined ranges of the NDVI values might be used for developing new fertilization strategies that are able to adjust N fertilization to wheat crop needs.


Introduction
One aim of agriculture has focused on the increase in cereal grain yields to feed an increasing population. Nitrogen (N) fertilization is an essential aspect for achieving this, as its application enhances grain yields. The optimization of the balance between yields and the environment needs fertilization strategies based on the regulation between N supply and crop needs [1]. N fertilizer should be applied to satisfy crop demand, taking into account the soil mineral N, the N mineralized from applied organic amendments (as farmyard manures, slurries, etc.), and the N mineralized from previous crop residues [2]. Besides, weather conditions vary considerably from growing season to growing season, releasing different amounts of N from organic matter, representing varying crop N needs within a single year throughout the season or among many years.
A key factor to obtain high yield is the optimal timing of N application [3]. The main Zadoks wheat growth stages [4] where topdressing N fertilizer is applied in the area where this experiment was performed (Araba, Basque Country, northern Spain) are the beginning of tillering (GS21) and the beginning of stem elongation (GS30). GS21 is the stage when tillers start emerging, where each one has the potential of producing a spike. Adequate N availability during the stem elongation stages determines the number of panicles and promotes crop growth [5,6]. In Araba, Basque Country, Spain, the mentioned growing stages are clues for establishing the required N fertilization, thus, usually, 40-60 kg N ha −1 is applied at GS21 and the greater and variable application is carried out at GS30. However, crop response to applied N is highly dependent upon the weather that occurs throughout the growing season [7]. The time elapsed between GS30 and grain harvest is long, and several factors might influence the N uptake by crops [8].
For achieving the maximum benefit with the minimum environmental impact, it is necessary to balance the maximum yield and the minimum cost for production. Regarding grain yield, some periods of N deficiency in cereals are critical and others have no impact [9]. Ravier et al. [10] showed that even intense N deficiencies in wheat crops may have no impact on grain yield as long as they occur around the early growing stages (GS30). However, N deficiencies should decrease in intensity later in the growing season to maintain satisfactory crop production. Jeuffroy and Bouchard [9] observed that the greatest losses in grain number were when the N deficiency was maintained until flowering (GS60 to GS69) [4]. Ravier et al. [10] detected threshold nitrogen nutrition index (NNI) values throughout key growing stages (GS30, GS32, GS39, and GS60), above which the yield did not decrease. However, NNI determination is time-consuming and needs destructive measurements of plant biomass and N content. Several authors have found correlations between NNI and optical sensing tool measurements [11][12][13]. In a previous paper, Aranguren et al. [14] studied if NNI could be predicted using RapidScan CS-45 or the chlorophyll meter Yara N-Tester TM , and they concluded that the NNI and the normalized difference vegetation index (NDVI) calculated from RapidScan CS-45) were well correlated in the humid Mediterranean conditions of Araba (R 2 = 0.70), opening the possibility of establishing threshold NDVI values throughout key wheat growing stages for achieving maximum yields. Using crop diagnostic tools is more practical than using NNI and offers real-time decision support tools for farmers. Besides, proximal sensing study might help in satellite-based remote sensing implantation. However, for accurate interpretation of satellite-based data necessary for making large-scale N fertilization recommendations, first, it is necessary to understand the information obtained in field-trials through proximal sensing, where different fertilization strategies are tested and indexes are measured at various growing stages, as done in the present study.
To achieve high yields without negative effects on the environment, it is necessary to synchronize N fertilization to wheat N needs. Thus, tracking the crop N status throughout the wheat growing season is key to finding the balance between crop N requirements and fertilizer needs. The present study aims to establish NDVI values for the soft winter wheat yield in humid Mediterranean conditions at the GS30, GS32, GS37, and GS65 growth stages with the RapidScan CS-45 instrument for the purpose of setting fertilizer requirements to maximize yields while avoiding N losses.

Study Site
Three field trials were established in Arkaute (Araba, Basque Country, Spain) at NEIKER facilities in three consecutive wheat growing seasons, defined as 2015, 2016, and 2017 in different fields under rainfed conditions (Figure 1). The three field trials were flat, presented similar characteristics, and the distances among them were 130 m. Soils had high pH values, were calcareous, presented a sandy clay loam texture, and had a moderate organic matter content (2-2.5%) in the first 0-30 cm. The soil was classified as Typic Calcexeroll [15]. Further details were described in a previous research study [7].

Climate
The climate of the area was humid Mediterranean according to the water regime of Papadakis' (1966) classification [16]. The meteorological conditions in the three growing seasons are described in Table 1.

Experimental Setup and Treatments
The experiment was of a factorial randomized complete block design with three factors (year, initial fertilization, and N rate at stem elongation) and featured four replicates. Three kinds

Measurements with RapidScan CS-45
In previous research, Aranguren et al. [14] found that NNI could be predicted using NDVI absolute values from the portable active crop canopy sensor RapidScan CS-45 (Holland Scientific, Lincoln, NE, USA). Following that finding, RapidScan CS-45 NDVI measurements from the key growing stages of the crop phenology (GS30, GS32, GS37, and GS65) were used with the aim of identifying and defining the NDVI threshold values that did not lead to a decrease in yield. Measurements were taken in all treatments in four replications and at control (0N) and overfertilized treatments (280N; Table 2).
The RapidScan CS-45 sensor features three optical photodetector channels that measure the reflectance from the crop at 670, 730, and 780 nm wavelengths. The RapidScan CS-45 sensors makes height not dependent on reflectance measurements, presenting a wide measurement range (from 0.3 m to 3 m). The field of view is approximately 45 degrees by 10 degrees. Measurements are sunlight independent, as RapidScan CS-45 is equipped with a modulated polychromatic lamp as an active light source. Besides, it is equipped with an internal GPS for geolocation with an accuracy < 1 m. It presents a low noise performance. The measurements with RapidScan CS-45 were taken as the sensor was passed over the crop surface at approximately 1 m, always maintaining a 90 • inclination with the crop at a constant walking speed. The sensor was handheld and two rows per plot were scanned. With both measurements, the coefficient of variation was calculated to measure the dispersion between readings. The coefficient of variation was always less than 5%. The NDVI values were averaged to generate a value from each plot.
The NDVI is the most widely known vegetation index (1) and is determined as follows: The accumulated NDVI ( NDVI) was calculated (2) for each treatment throughout the wheat growing season with the measurements taken at GS30, GS32, GS37, and GS65 as follows: Remote Sens. 2020, 12, 2749 5 of 17

Statistical Analysis
Two different statistical analyses were carried out: different analyses of variance (ANOVAs) were performed for the factors affecting wheat grain yield and NDVI, and a coefficient of determination (R 2 ) was performed for the relationship between NDVI and wheat grain yield. The yield-influencing three factors analyzed were the growing season, initial fertilization treatment, and N rate at GS30, determined via analyses of variance (ANOVA) using the R 3.2.5 software package. There was an interaction among the growing season, initial fertilization treatment, and N rate at GS30. Therefore, it was necessary to analyze each factor depending on the other factors. An ANOVA was performed to analyze differences among the N rate at GS30 in each initial fertilization treatment and each growing season (2015, 2016, and 2017). Another ANOVA was performed to analyze the differences among the initial fertilization treatment in each N rate at GS30 and growing season (2015, 2016, and 2017). To separate the means, Duncan's test was used (p < 0.05), using the R package "agricolae" [17]. 0N and 280N were not included in the ANOVA. The statistical significance of the difference in yield between each treatment and 0N, as well as each treatment and 280N, were determined according to Welsh's t-test (for unequal variances) or via a pooled variance t-test (for equal variances). Levene's test was carried out with the aim of determining if the variance in the two populations to be tested was equal or not. The results of the t-tests are presented in Table A1 (differences with 0N) and Table A2 (differences with 180N).
As in the case of yield, it was necessary to analyze each factor depending on the others. The same was done for analyzing differences in the NDVI, since the aim was to relate NDVI values with yield. An ANOVA was performed to analyze the differences in the NDVI values among the N rate applied at GS30 in each initial fertilization treatment and growing season (2015, 2016, and 2017) using R 3.2.5. Another ANOVA was performed to analyze the differences among initial fertilization treatment in each N rate at GS30 and growing season (2015, 2016, and 2017). To separate the means, Duncan's test was used (p < 0.05) using the R package "agricolae" [17]. The letters from Duncan's tests were used to see the differences among means and are not presented in the tables (Tables 3-5) because it would make their visualization difficult and heavy. Here, letters are presented in Tables A3 and A4. Nevertheless, when differences among NDVI values of N fertilization treatments were significant, each growing stage was marked with *** (i.e., significant at a 0.001 probability level) in Tables 3-5. When there was no difference, each growing stage was marked with "ns." Here, 0N and 280N were not included in the ANOVA. As with the yield, the statistical significance of the differences in NDVI between each treatment and 0N, as well as each treatment and 280N, were determined according to Welsh's t-test (unequal variances) or a pooled variance t-test (equal variances). Levene's test was carried out with the aim of determining if the variance in the two populations to be tested was equal or not. The results obtained from the t-tests are presented in Table A1 (differences with 0N) and Table A2 (differences with 180N).
A coefficient of determination (R 2 ) for the relationship between NDVI and wheat grain yield was performed using R 3.2.5.

Period Comprehended from GS21 to GS30
The NDVI values at GS30 in 2016 (Table 4) were higher than the ones in 2015 and 2017 (Tables 3  and 5), especially for treatments with organics as initial fertilizers. In 2015 and 2017, the values at GS30 were different for three initial fertilization treatments, where the values at the conventional treatment were higher (0.55-0.60) than those of the organic treatments (0.40-0.45). However, in 2016, the values for three initial fertilization treatments at GS30 were similar (0.65). In 2016, even the control (0N) reached a value higher than 0.5. More time elapsed in 2016 from GS21 to GS30 than in 2015 and 2017 (Table 2). Besides, in 2016 it rained far less before GS21 (Table 2) than in 2015 and 2017. Regarding the overfertilized plots, the NDVI value in 2015 was 0.6, in 2016 was 0.75, and in 2017 was 0.7.

Period Comprehended from GS30 to GS32
In 2016, for the treatments where 0 kg N ha −1 was added at GS30, the values descended from GS30 to GS32 (Table 4). In 2015, following the nitrogen application at GS30, the NDVI values rapidly increased, showing an increase in the photosynthetically active biomass. Besides, the 0N values also increased but in a smaller degree. In 2015 (Table 3) at GS32 in conventional treatment, there were no differences in the NDVI values when N was applied at GS30 (Table A3). In treatments where organics were applied, there were no differences among the highest N rates (80, 120, 160 kg N ha −1 ). In 2016 (Table 4), at GS32, more differences among treatments were appreciable according to the N rate applied at GS30 (the higher the N rate, the higher the NDVI value). In 2017, the NDVI values did not increase from GS30 to GS32 in any of the treatments and there were no differences among the N rates applied at GS30 (Table 5). Regarding the NDVI values in the overfertilized plot in 2017, they decreased from GS30 (0.7) to GS32 (0.65).

Period Comprehended from GS32 to GS37
In 2015, the NDVI values showed that from GS32 onwards the organic treatments and conventional treatment with 80, 120, 160 kg N ha −1 rates applied at GS30 were similar (Table 3), reaching values above 0.75. In 2016, the NDVI values differed according to the different rates applied at GS30 in three initial fertilization treatments from GS32 on ( Table 4). In 2017, the different N rates applied at GS30 were not reflected by different NDVI values in any of the initial fertilizer treatments. However, differences in the NDVI values between conventional treatment and organic treatments persisted until GS37 (Table 5). In 2017, in the overfertilized plot, the NDVI values remained much higher than the rest of treatments (0.7).

Period Comprehended from GS37 to GS65
In 2015, the NDVI values in the majority of treatments decreased from GS37 to GS65 (Table 3). In 2016, the NDVI values in all treatments remained the same from GS37 to GS65 (Table 4). In 2017, the differences among NDVI depending on the N rates applied at GS30 were not measurable until GS65 in three initial fertilization treatments (Table 5). Thus, the higher the N rate applied at GS30, the higher the NDVI values at GS65. The NDVI values at the conventional treatments remained higher than the organics for the lower N rates applied at GS30 (0 and 40 kg N ha −1 ). Table 3. Wheat grain yield (kg N ha −1 ) and normalized difference vegetation index (NDVI) readings collected with RapidScan CS-45 in 2015 at stem elongation, second node, leaf-flag emergence, and mid-flowering (GS30, GS32, GS37, and GS65, respectively (Zadoks et al. 1974)) in the field experiment located in Arkaute.  Table 3. Wheat grain yield (kg N ha −1 ) and normalized difference vegetation index (NDVI) readings collected with RapidScan CS-45 in 2015 at stem elongation, second node, leaf-flag emergence, and mid-flowering (GS30, GS32, GS37, and GS65, respectively (Zadoks et al. 1974)) in the field experiment located in Arkaute.  Overfert.

280N nd nd
ANOVAs (analyses of variance) were performed to analyze differences in the yield and NDVI values (at GS32, GS37, and GS65) among N rates applied at GS30 (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) in each initial fertilization treatment (conventional, dairy slurry, and sheep manure). Different capital letters represent differences among yields at different N application treatments at GS30 for each initial fertilization treatment. Differences from Duncan's test in NDVI values are shown in Table A3 and Table A4. ***, significant at a 0.001 probability level. 0N was not included in the ANOVAs. The statistical significance of the differences in yield and NDVI between each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 0N are shown in Table A1 and  Table A2. sd, standard deviation. nd, no data. Overfert.

280N nd nd
ANOVAs (analyses of variance) were performed to analyze differences in the yield and NDVI values (at GS32, GS37, and GS65) among N rates applied at GS30 (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) in each initial fertilization treatment (conventional, dairy slurry, and sheep manure). Different capital letters represent differences among yields at different N application treatments at GS30 for each initial fertilization treatment. Differences from Duncan's test in NDVI values are shown in Table A3 and Table A4. ***, significant at a 0.001 probability level. 0N was not included in the ANOVAs. The statistical significance of the differences in yield and NDVI between each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 0N are shown in Table A1 and  Table A2. sd, standard deviation. nd, no data. ANOVAs (analyses of variance) were performed to analyze differences in yield and NDVI values (at GS32, GS37 and GS65) among N rates applied at GS30 (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) in each initial fertilization treatment (conventional, dairy slurry and sheep manure). Different capital letters represent differences among yields at different N application treatments at GS30 for each initial fertilization treatment. Different lowercase letters represent differences among yields for different initial fertilization treatment at the same N application treatment at GS30. Differences from Duncan's test in NDVI values are shown in Table A3 and Table A4. ns, not significant. ***, significant at a 0.001 probability level. Here, 0N and 280N were not included in the ANOVAs. The statistical significance of the differences in yield and NDVI between each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 0N, and each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 280N are shown in Table A1 and Table A2. sd, standard deviation.
ANOVAs (analyses of variance) were performed to analyze differences in the yield and NDVI values (at GS32, GS37, and GS65) among N rates applied at GS30 (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) in each initial fertilization treatment (conventional, dairy slurry, and sheep manure). Different capital letters represent differences among yields at different N application treatments at GS30 for each initial fertilization treatment. Differences from Duncan's test in NDVI values are shown in Tables A3 and A4. ***, significant at a 0.001 probability level. 0N was not included in the ANOVAs. The statistical significance of the differences in yield and NDVI between each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 0N are shown in Tables A1 and A2. sd, standard deviation. nd, no data.   ANOVAs (analyses of variance) were performed to analyze differences in yield and NDVI values (at GS32, GS37 and GS65) among N rates applied at GS30 (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) in each initial fertilization treatment (conventional, dairy slurry and sheep manure). Different capital letters represent differences among yields at different N application treatments at GS30 for each initial fertilization treatment. Differences from Duncan's test in NDVI values are shown in Table A3 and Table A4. ***, significant at a 0.001 probability level. Here, 0N and 280N were not included in the ANOVAs. The statistical significance of the differences in yield and NDVI between each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 0N, and each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 280N are shown in Table A1 and Table A2  ANOVAs (analyses of variance) were performed to analyze differences in yield and NDVI values (at GS32, GS37 and GS65) among N rates applied at GS30 (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) in each initial fertilization treatment (conventional, dairy slurry and sheep manure). Different capital letters represent differences among yields at different N application treatments at GS30 for each initial fertilization treatment. Differences from Duncan's test in NDVI values are shown in Table A3 and Table A4. ***, significant at a 0.001 probability level. Here, 0N and 280N were not included in the ANOVAs. The statistical significance of the differences in yield and NDVI between each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 0N, and each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 280N are shown in Table A1 and Table A2  ANOVAs (analyses of variance) were performed to analyze differences in yield and NDVI values (at GS32, GS37 and GS65) among N rates applied at GS30 (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) in each initial fertilization treatment (conventional, dairy slurry and sheep manure). Different capital letters represent differences among yields at different N application treatments at GS30 for each initial fertilization treatment. Different lowercase letters represent differences among yields for different initial fertilization treatment at the same N application treatment at GS30. Differences from Duncan's test in NDVI values are shown in Table A3 and Table A4. ns, not significant. ***, significant at a 0.001 probability level. Here, 0N and 280N were not included in the ANOVAs. The statistical significance of the differences in yield and NDVI between each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 0N, and each treatment (X + 0N, X + 40N, X + 80N, X + 120N, ANOVAs (analyses of variance) were performed to analyze differences in yield and NDVI values (at GS32, GS37 and GS65) among N rates applied at GS30 (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) in each initial fertilization  Tables A3 and A4. ***, significant at a 0.001 probability level. Here, 0N and 280N were not included in the ANOVAs. The statistical significance of the differences in yield and NDVI between each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 0N, and each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 280N are shown in Tables A1 and A2. sd, standard deviation. Table 5. Wheat grain yield (kg N ha −1 ) and NDVI readings collected with RapidScan CS-45 in 2017 at stem elongation, second node, leaf-flag emergence and mid-flowering (GS30, GS32, GS37, and GS65, respectively (Zadoks et al. 1974)) in the field experiment located in Arkaute.  ANOVAs (analyses of variance) were performed to analyze differences in yield and NDVI values (at GS32, GS37 and GS65) among N rates applied at GS30 (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) in each initial fertilization treatment (conventional, dairy slurry and sheep manure). Different capital letters represent differences among yields at different N application treatments at GS30 for each initial fertilization treatment. Different lowercase letters represent differences among yields for different initial fertilization treatment at the same N application treatment at GS30. Differences from Duncan's test in NDVI values are shown in Table A3 and Table A4. ns, not significant. ***, significant at a 0.001 probability level. Here, 0N and 280N were not included in the ANOVAs. The statistical significance of the differences in yield and NDVI between each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 0N, and each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 280N are shown in Table A1 and Table A2  ANOVAs (analyses of variance) were performed to analyze differences in yield and NDVI values (at GS32, GS37 and GS65) among N rates applied at GS30 (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) in each initial fertilization treatment (conventional, dairy slurry and sheep manure). Different capital letters represent differences among yields at different N application treatments at GS30 for each initial fertilization treatment. Different lowercase letters represent differences among yields for different initial fertilization treatment at the same N application treatment at GS30. Differences from Duncan's test in NDVI values are shown in Table A3 and Table A4. ns, not significant. ***, significant at a 0.001 probability level. Here, 0N and 280N were not included in the ANOVAs. The statistical significance of the differences in yield and NDVI between each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 0N, and each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 280N are shown in Table A1 and Table A2  ANOVAs (analyses of variance) were performed to analyze differences in yield and NDVI values (at GS32, GS37 and GS65) among N rates applied at GS30 (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) in each initial fertilization treatment (conventional, dairy slurry and sheep manure). Different capital letters represent differences among yields at different N application treatments at GS30 for each initial fertilization treatment. Different lowercase letters represent differences among yields for different initial fertilization treatment at the same N application treatment at GS30. Differences from Duncan's test in NDVI values are shown in Table A3 and Table A4. ns, not significant. ***, significant at a 0.001 probability level. Here, 0N and 280N were not included in the ANOVAs. The statistical significance of the differences in yield and NDVI between each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 0N, and each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 280N are shown in Table A1 and Table A2  ANOVAs (analyses of variance) were performed to analyze differences in yield and NDVI values (at GS32, GS37 and GS65) among N rates applied at GS30 (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) in each initial fertilization treatment (conventional, dairy slurry and sheep manure). Different capital letters represent differences among yields at different N application treatments at GS30 for each initial fertilization treatment. Different lowercase letters represent differences among yields for different initial fertilization treatment at the same N application treatment at GS30. Differences from Duncan's test in NDVI values are shown in Table A3 and Table A4. ns, not significant. ***, significant at a 0.001 probability level. Here, 0N and 280N were not included in the ANOVAs. The statistical significance of the differences in yield and NDVI between each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 0N, and each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 280N are shown in Table A1 and Table A2. sd, standard deviation.
ANOVAs (analyses of variance) were performed to analyze differences in yield and NDVI values (at GS32, GS37 and GS65) among N rates applied at GS30 (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) in each initial fertilization treatment (conventional, dairy slurry and sheep manure). Different capital letters represent differences among yields at different N application treatments at GS30 for each initial fertilization treatment. Different lowercase letters represent differences among yields for different initial fertilization treatment at the same N application treatment at GS30. Differences from Duncan's test in NDVI values are shown in Tables A3 and A4. ns, not significant. ***, significant at a 0.001 probability level. Here, 0N and 280N were not included in the ANOVAs. The statistical significance of the differences in yield and NDVI between each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 0N, and each treatment (X + 0N, X + 40N, X + 80N, X + 120N, X + 160N) and 280N are shown in Tables A1 and A2. sd, standard deviation.

NDVI Values for Maximum Grain Yield
Grain yield was significantly influenced by the amount of mineral N fertilizer at GS30 in the three growing seasons. Different NDVI combinations for similar yields were detected throughout the three growing seasons (Tables 3-5). It is remarkable that yields in 2016 were higher than in the other two growing seasons, and in 2017 yields were especially low (with the exception of the overfertilized plot (8000 kg ha −1 )). Overall, and taking into account the controls, wheat grain yield varied between 4100 and 8700 kg ha −1 in 2015 (Table 3), 5200-10,700 kg ha −1 in 2016 (Table 4) and 3350-8000 kg ha −1 in 2017 (Table 5). In 2015, the maximum wheat yield was achieved in all the initial fertilization treatments (conventional, dairy slurry, and sheep manure) with 80 kg N ha −1 at GS30. Regarding the NDVI values, the treatments with 80, 120, and 160 kg N ha −1 rates applied at GS30 were comparable (Tables 3 and A3) in three initial fertilization treatments reaching values above 0.75, 0.76, and 0.69 at GS32, GS37, and GS65, respectively. The absence of significant differences in the NDVI values among treatments with the 80, 120, and 160 kg N ha −1 rates matched with the N rate with which the maximum yield was achieved (80 kg N ha −1 at GS30). In fact, for treatments with lower N application rates at GS30 (0 y 40 kg N ha −1 ), the NDVI values (NDVI = 0.50-0.70) did not reach the 2015 NDVI values mentioned for the high yield producing treatments. In 2016 and 2017 (Tables 4 and 5), with conventional and dairy slurry treatments, the maximum wheat yield was achieved with a 80 kg N ha −1 rate and in sheep manure treatment with a 120 kg N ha −1 rate at GS30. In 2016, the NDVI values at the treatments that achieved the highest yield were 0.70 at GS32 and 0.75 at GS37 and GS65 in three initial fertilization treatments (Table 4). In the treatments where the highest yields were not achieved, the NDVI values did not reach the mentioned 2016 values (0.70 at GS32 and 0.75 at GS37 and GS65) corresponding to the most productive treatments. In 2017, the organics applied as initial fertilizers produced less than the conventional treatments for each N rate applied at GS30. In conventional treatment, the NDVI values at GS32 and GS37 were 0.5 for all different N rates applied at GS30. In the case of organics as initial fertilizers, the NDVI values at GS32 and GS37 were around 0.4 for all different N rates applied at GS30. The higher NDVI values achieved in conventional treatments than in organic treatments at those growing stages (Table A4) meant a difference of 500-1500 kg ha −1 in yield, depending on the N rate applied at GS30. At GS65, the NDVI values at the treatments that achieved the highest yield in three types of initial fertilization were 0.65 (Table 5), while in the rest of treatments the NDVI remained lower (0.4-0.6). That increase at GS65 meant a difference of up to 2000 kg ha −1 in yield between the treatments with the maximum yields and the lowest yields (the ones that received 0 kg N ha −1 at GS30). In 2017, in the overfertilized plot, the NDVI values remained higher than the rest of treatments from GS30 to GS65 (0.7 at GS30, 0.65 at GS32, 0.7 at GS37, and 0.7 at GS65; Table A2), which meant a difference in yield of 1500 kg ha −1 in comparison with maximum yields in conventional treatment, and 2000-2500 kg ha −1 in comparison with the maximum yields in the organic treatments.
Different NDVI combinations for similar yields were detected throughout the three growing seasons (Tables 3-5

Discussion
These results show the high potential of the RapidScan CS-45 instrument to track the N status of wheat plants with NDVI during the crop growth, as Li et al. [13] have also shown for rice. In fact, the NDVI has been considered a good indicator of green biomass or N content, focusing on the plant canopy [18], and has been a widely used vegetation index in agriculture for the last 40 years

Discussion
These results show the high potential of the RapidScan CS-45 instrument to track the N status of wheat plants with NDVI during the crop growth, as Li et al. [13] have also shown for rice. In fact, the NDVI has been considered a good indicator of green biomass or N content, focusing on the plant canopy [18], and has been a widely used vegetation index in agriculture for the last 40 years [19]. For following the physiological processes that control plant N uptake and yield, real-time information about crop nutritional status is necessary [6]. Each phenological period has an influence on the wheat harvest metrics, and thus the N deficiencies present a different effect depending on the intensity and duration in each moment of the crop growth [9]. Some periods of N deficiency in wheat are critical and others have no impact or are tolerable [10], making yields variable among years. In the present study, different NDVI evolution patterns were detected, causing different effects on grain yield (Tables 3-5). Ravier et al. [10] detected threshold NNI values throughout key growing stages that were tolerable for wheat without decreasing yields. NNI values should always exceed 0.4, 0.7, 0.7, and 0.8 at growth stages GS30, GS32, GS39, and GS60, respectively.

NDVI Values at Key Growing Stages
In the present study the minimum NDVI values for achieving the highest yields were 0.7-0.75 at GS32, GS37, and GS65. However, the NDVI threshold values at GS30 were not clear (0.4-0.65). In a previous paper [14], good correlations between the NNI and NDVI were found; therefore, those correlations were used to compare the results of this study with the ones published by Ravier et al. [10]. For achieving the NNI values proposed by Ravier et al. [10], the NDVI values should be 0.5 at GS30, 0.75 at GS32, 0.75 at GS37, and 0.8 at GS65. Similarities and differences were found when comparing the results of the present study with the ones proposed by Ravier et al. [10]. At GS30, the results obtained in the present study and the ones proposed by Ravier et al. [10] diverged. However, at GS32, GS37, and GS65, the values obtained from this study were similar to the ones obtained by Ravier et al. [10].
Focusing on the GS30 growing stage, in all treatments in 2016 and only for the conventional treatments in 2017 and 2015 were the NDVI values higher than the ones proposed by Ravier et al. [10], 0.4 for the NNI and 0.5 for the NDVI according to the relationship established by Aranguren et al. [14]. The same happened for all overfertilized treatments, where the NDVI values were 0.75, 0.70, and 0.6 in 2016, 2017, and 2015, respectively. Tillering (GS21 to GS29) has been reported to be important for reaching an adequate wheat yield [6,20], as a good tiller establishment allows a greater potential for biomass accumulation [6]. At this early growing stage (before GS30), when the soil water content is high, the most limiting contributor to slow tillering could be the soil temperature and anaerobic conditions [6]. The tillering period is related with wheat yield as it is involved in the determination of the spike number [21]. Prolonging the tillering period results in a higher number of tillers formed [6], which is probably what happened in 2016. The time elapsed from GS21 to GS30 was longer in 2016 (56 days) than in 2017 (35 days) and in 2015 (30 days) (Table 1), promoting more efficient N incorporation (higher N uptake), leading to a higher photosynthetically active biomass. The longer the time elapsed from GS21 to GS30, the higher the achieved NDVI values (0.65 in 2016, 0.45-0.6 in 2017, and 0.38-0.55 in 2015). Otherwise, in 2015 and 2017, it rained much more than in 2016, with a difference of 353 mm and 103 mm, respectively, from the same period in 2016 (Table 2). These heavy rains caused soil waterlogging, visible puddles, and anaerobic conditions, explaining the lower NDVI values. Soil waterlogging in winter in this area is not uncommon in very wet years. In fact, the experimental area is located in the Quaternary aquifer of Vitoria-Gasteiz, where the existence of several historically built trenches that 0.5 m deep with the purpose of avoiding flooding in agricultural areas have been described [22]. Besides, those environmental conditions have provoked less N availability from the organics applied as initial fertilizers, explaining the lower NDVI values in the treatments with organics in 2015 and 2017 (0.4-0.45) than the conventional treatments (Tables 3 and 5). Otherwise, those values were lower than the threshold values proposed by Ravier et al. [10], 0.4 for NNI and 0.5 for NDVI according to the relationship established by Aranguren et al. [14]. In 2015, the yields were comparable between organic treatments and conventional treatments for each N rate applied at GS30. That fact suggests that a N deficiency at GS30 can be partially corrected if the plant is able to absorb enough N in the following growing stages, as in the year 2015 (0.4, 0.7, 0.7 (organics) and 0.55, 0.7, 0.7 (conventional) for GS30, GS32, and GS37, respectively). Similarly, Morris et al. [23] reported that a N deficiency in early growing stages results in maximum or near maximum yields in winter wheat. Conversely, in 2017, the yields in organic treatments were significantly lower than the yields in conventional treatments for the same N rate applied at GS30, as the NDVI values in conventional treatments remained higher than in organics, at least until GS37 (Table A4; discussed in the following paragraphs). Yields in 2016 were generally higher than in 2015 and 2017, suggesting that the period before GS30 always has an important effect on crop yield if it maintained in the following key growing stages (from GS32 to GS65). No specific value was found for the GS30 growing stage, as the data volume was small for all possible options presented in the following key wheat growing stages (GS32, GS37, and GS65). That high variability made the establishment of the GS30 value difficult, thus, new data obtained from more experiments might help in its adjustment. Ravier et al. [10] also remarked that the N threshold at GS30 was not so clear in terms of the comparison with the other stages.
Focusing on stem extension (GS32 to GS37), the values detected in the present study for achieving maximum yields and the ones proposed by Ravier et al. [10] were similar (a NDVI value around 0.75). Adequate N availability at stem extension promotes crop growth through the development of viable tillers, increasing the future sink capacity of the plant, being related to the spike numbers and the grain number per unit area [9]. A N deficiency in this period leads to the greatest losses in grain number [9]. It is relevant that in 2017, during this period, the NDVI values were very low ( Table 5) in most of the treatments (0.4-0.55) because of the lack of rain after N application at GS30 (Table 1), justifying the lower yields achieved in 2017 in comparison with 2015 and 2016. In fact, the treatments with organics presented lower NDVI values (NDVI = 0.40-0.45) than the conventional treatments (NDVI = 0.52-0.55), explaining the significantly higher yields in conventional treatments, as mentioned in the previous section. The NDVI values were significantly higher in the 280N treatment (Table A2, NDVI = 0.65-0.70) because of the high N rate applied at GS21 (80 kg N ha −1 ), explaining the higher yields achieved (8020 kg ha −1 ). Anyway, these values were lower than the threshold values proposed by Ravier et al. [10] and the ones detected in the present study (NDVI around 0.75). Villegas et al. [24] showed that drought stress at stem extension resulted in lower rates of biomass accumulation. Adequate water availability for crops after the application of N fertilizer can be beneficial for a good N status and for achieving high yields [25], as in 2015 and 2016 (Table 1), where after the mineral fertilizer application at GS30 the soil was sufficiently moist, making the N uptake possible. The rainfall quantity necessary for the absorption of a N application by the crops is at least 15 mm in the fortnight after the application [26].
At the GS65 growing stage, the NDVI values achieved at the treatments with the highest yields were similar to the ones proposed by Ravier et al. [10], with NDVI values around 0.75. In the 2017 growing season, this was the first growing stage where differences in the NDVI values among different N rates were appreciable ( Table 5). The previous dry period and the rainfall that occurred after GS37 made the uptake of the N applied at GS30 viable. The N absorption by the crop from GS37 to GS65 was evident when comparing the yields obtained at the lower and higher doses applied at GS30, although it has been reported by other authors that this is late for increasing yields [27]. Similar to the present study, Asseng et al. [28] reported that it is a viable method to improve yields via nitrogen application at GS59 (before flowering) if the soil is wet. Anyway, in 2017, there was a smaller increase in yield per kg of N applied at GS30 in comparison with the 2015 and 2016 growing seasons, which is appreciable in the lower yields achieved and in the narrow difference of up to 2000 kg ha −1 in yield between the treatments with the maximum yields and the lowest yields.

NDVI Dynamics and Wheat Grain Yield
In the present study, it was shown that when values at GS30 were 0.55-0.6 (higher than the ones proposed by Ravier et al. [10]) and around 0.75 at GS32, GS37, and GS65, maximum yields were achieved. It is remarkable that different dynamics in wheat NDVI values throughout the growing cycle were possible for achieving a specific yield, as mentioned in the results. There were treatments with a low N status at GS30 that raised their NDVI values by GS32 or GS37, while there were other treatments that, due to the low N rate applied at GS30, reduced their NDVI values by GS32 or GS37. During the wheat growing cycle, the circumstances in the early, middle, and late growth stages are strongly associated with the number of tillers, grain per spike, and 1000 grain weight, respectively [29]. A low N status during a given growing period, thus causing a deficiency in one of these parameters, may not have an effect on yield if the N status is higher in another key period. Thus, a nitrogen deficiency that may affect the number of tillers can be supplied by a higher N status when the grains per spike are established [9]. The use of NDVI is a way of testing that every key growing stage (GS30, GS32, GS37, and GS65) has an important effect on crop yield. The NDVI from early growing stages to anthesis holds a large amount of crop growth data, thus representing the impact of plant growth status on yield formation in wheat, as Wang et al. [29] have shown. Besides, the NDVI from GS30 to GS65 better explained the grain yield variability (80% of the variability, Figure 2) than just using the measurement of one growing stage alone. In a previous paper, Aranguren et al. [14] found that NDVI absolute values could explain yields at GS65 (70% of the variability). In the case that a period of dryness occurs after the N application, as in 2017, where the N absorption did not occur until late growing stages, by using the NDVI it is possible to predict how this will affect the wheat grain yield. Saturation effects have been reported for NDVI values around 0.8 when NNI > 0.8 [14]. Anyway, the NNI threshold values for achieving maximum yields in wheat proposed by Ravier et al. [10] were always a NNI value of 0.8 or lower. Therefore, in this case, the saturation effect of the NDVI would not be problematic. Cao et al. [30] reported that saturation effects can be decreased by using NDRE (a red-edge-based vegetation index), whereas in studies with the RapidScan CS-45 instrument it was shown that the NDRE reaches the same precision as the NDVI and a comparable saturation level [14,31]. Besides, thinking of a possible future application of the images from the Sentinel satellite, it should be considered that the NDRE has a lower resolution than NDVI, as the red edge band (used for NDRE calculation) and red band (used for NDVI calculation) have radiometric resolutions of 20 m and 10 m, respectively [32].

Proximal Sensing for Correcting Wheat Yield
For designing N fertilization schedules related to specific agronomic regions, it is necessary to obtain better knowledge of the critical moments and critical values for adjusting the requirements of the N fertilizer by the crop during the key vegetative growing stages. Remote sensing measurements are not laborious and can be taken periodically to monitor crop N status in an effective way until the end of the cereal-growing season. The present results show that following the NDVI dynamics throughout crop growth would be helpful for determining if there is a N deficiency, if it is possible to correct it by applying N, and therefore optimize crop N management for achieving high yields. In fact, it has been reported that it is always possible to correct a nitrogen deficiency until flowering [10] if soil is humid [33]. These NDVI threshold values have been defined for indicated edaphoclimatic conditions here and for the Cezanne variety. The fact that the results reported in the present study are similar to those obtained by Ravier et al. [10] in France with 215 treatments, different wheat varieties, and different locations, suggest that these values might be used in a wider range of conditions. Caution is needed regarding the universality of minimum NDVI values across wheat varieties and geographical locations. As Diacono et al. [1] have reported, wheat variety should be considered, as each variety might have its own calibration requirements. They recommended the use of data normalization, but several concerns related to the use of normalized data have been found. Aranguren et al. [14] found that the relationship between the NNI and normalized NDVI values is worse than using absolute NDVI values. Ravier et al. [12] reported that it is not easy to find a control strip to be representative of the entire field and ensure that an overfertilized fringe is not N deficient. Besides, data normalization makes the use of proxy tools more complicated for farmers. Bonfil et al. [31] reported that using the RapidScan CS-45 instrument for cultivar classification, combining the data from a few monitoring days, was possible.
Understanding ground-based remote sensing could help with satellite remote sensing implementation, and this is why accurate interpretations of satellite-based data requires robust ground-based reference data [34]. However, satellite remote sensing requires special training to process data, while active canopy sensors are simple to operate but are not adequate for extended surfaces.

Conclusions
The threshold NDVI values measured with a RapidScan CS-45 instrument for achieving the highest yields in the wheat Cezanne variety under humid Mediterranean conditions were 0.7-0.75 at the GS32, GS37, and GS65 growing stages. It was not clearly established here for the GS30 growing stage, as N deficiencies (NDVI = 0.4-0.45) may not affect wheat yields as long as the N status increases by GS32 and it is maintained thereafter.
Following the NDVI dynamics throughout the growing season could help to predict the yields at harvest time. Here, the ΣNDVI from GS30 to GS65 explained 80% of wheat yield variability. Therefore, a given yield could be achieved with different dynamics in wheat NDVI values throughout the growing cycle. The established ranges of NVDI values might be used for adjusting N fertilization to the wheat crop needs.