Remote and Proximal Sensing Applications for Durum Wheat Nutritional Status Detection in Mediterranean Area

Combining remote and proximal sensing in agriculture is essential to monitor crop spatialtemporal variability and to provide high-quality prescription maps for the precision agriculture applications. The study showed how different combinations of soil management (no tillage—NT vs. conventional tillage—CT) and nitrogen (N) fertilization levels (0.90 and 180 kg N ha−1) can affect the durum wheat nutritional status and development through vegetation indices computation and proximal sensing tool application. Chlorophyll and N crop content were measured, in addition a proximal sensing tool and multispectral imagery equipped on unmanned aerial vehicle were used. The N input is the key driver for durum wheat development (4.5 ± 0.92 t ha−1 on average), but when it was not provided the NT performed better than CT (2.51 ± 0.22 vs. 1.46 ± 0.28 t ha−1 respectively) in terms of grain yield. This is due to the greater content of organic matter and N availability which characterizes the NT system. The near infrared (NIR) band-based vegetation indices can well detect the durum wheat nutritional status (R2 = 0.70 on average). The showed results can provide an important contribution in the implementation of ago-environmental policies aimed at environmental impact of cereal-based-cropping systems reduction.


Introduction
In the European Union, the production of wheat is estimated at 152 million kg 10 3 , about 10% above the harvest of 2018 [1]; in detail in Mediterranean semi-arid climate regions durum wheat (Triticum turgidum L. subsp. durum) is among the most common winter cereal crops [2,3]. In this geographic area drought stress as a combination of erratic rainfall distribution and high temperature during grain filling stages, is the main constraint limiting crop production [4] and determining their instability.
In the near future, the Mediterranean area is expected to face more severe drought and an increase in average temperature, due to climate change [5]; climate change and the simultaneous global population increase, estimated at 2.6 billion by 2050 [6], require a range of environmental measures such as optimizing resources [7], but production must remain high and environmental degradation must be reduced at the same time. If modern agriculture with conventional tillage techniques have allowed high yields, it is also equally true that such soil management induced greater soil erosion [3].
Both the reduced soil fertility and the effects of climatic change strongly affect durum wheat yields and quality traits [8]; to counteract such criticality, conservative agriculture (CA) practices-based on minimum or no soil disturbance, permanent residue cover and planned crop rotations or associations [9] should be applied; these aim to promote water conservation, biological processes, and organic matter accumulation due to reduced soil erosion and CO 2 emissions [10][11][12][13].
However, not all studies on CA practices have led to the same results regarding the crop yield. In a dry context, the CA can have a positive effect on crop productivity due to increased soil moisture and nutrient availability [14]. In the literature, there are often contradictory results when the long-term impact of CA practices for durum wheat yield is analysed.
This is due to differences among experimental sites, such as soil and weather conditions, different experimental designs, and different agronomic management of the site [15,16] for example [17] reported 20 years of experiment results on durum wheat grain yield, and concluded that for durum wheat, NT did not allow substantial yield benefits leading to comparable yields with respect to CT in 10 out of 20 years.
However, these studies have always focused on the final grain yield and not on the crop cycle or nutritional status. This can be easily monitored using a precision agriculture (PA) instrument like remote sensing.
In fact, PA allows the use of farm inputs such as fertilizers and herbicides to be optimized [18,19], and thus maximize profit and minimize negative environmental impacts [16] by addressing spatial variability.
PA methods hold great potential as a tool for: high-throughput phenotyping for plant breeding [20,21], decision making for precision agriculture [22,23], predicting yields [24] and predicting spatial field variability in experimental sites [25].
In order to apply PA methods, spatial-temporal crop variability must be measured, and to do so remote and proximal sensing tools can be used and combined [26]. Their usefulness relies on the fact that they are non-destructive, low time-consuming methods, moreover they are also well-correlated with agronomical and important physiological crop traits [27].
An example of those sensing tools is the unmanned aerial vehicle (UAV) equipped with multispectral camera. This allows to monitor the crop with high spatial-temporal resolution during the production season [28].
Spectral remote sensing is a widely used tool for agricultural production and crop management with respect to site-specific fertilizer applications [29,30]; its potential as a vehicle-based high-throughput phenotyping technology under field conditions has been demonstrated [31][32][33][34].
Studies regarding the effect of CA practices and several nitrogen fertilization levels on chlorophyll content and nitrogen status of durum wheat in relation to vegetation indices are lacking. For example, some studies showed that the yield estimated from the normalized difference vegetation index (NDVI) had a strong relationship with grain yield in wheat [35,36] in addition, NDVI has also been used to estimate crop growth status based on the different patterns of reflection of green organs and soil in wheat and other cereals [37][38][39][40].
Following the approach of [41], it is necessary implement the knowledge of the relationships between field data (such as chlorophyll and nitrogen content in the crop, yield components as the number of kernels per spike) and remote-sensing data (such as several vegetation indices) in relation to the type of management operated on durum wheat.
Based on the above considerations, the aims of this study are to evaluate how different soil management combinations and nitrogen fertilization levels can affect: • the durum wheat nutritional status, chlorophyll content and yield components; • the effectiveness of proximal and remote sensing tools to well correlate the durum wheat chlorophyll content and nutritional status; in addition, the scientific findings emerging from this experiment will be used to provide a contribution in the implementation of the agri-environmental policies adopted in the examined context.

Experimental Site
The experimental site is located in the "Pasquale Rosati" experimental farm of the Polytechnic University of Marche in Agugliano, Italy (43 • 32 N, 13 • 22 E, at an altitude of 100 m above sea level and a slope gradient of 10%) ( Figure 1) on a Silty-clay soil classified as Calcaric Gleyic Cambisol [42]. The weather data were provided by ASSAM (Agenzia Servizi al Settore Agroalimentare delle Marche, Osimo stazione, Ancona, Italy). The thermo-pluviometry trend data (Table 1) showed a similar trend in both temperatures and rainfall with a slight difference of 0.4 • C in the mean air temperature and 30 mm of rainfall between the years 2018 and 2019.  2018  29  173  143  37  95  48  57  25  19  67  42  61  796  2019  70  22  36  59  165  1  105  15  112  62  66  53  765 Soil sampling was performed with a hand auger (model: Suelo HA-3, Zhejiang Lujian Instrument Equipment Co., Ltd., Zhejiang, China) immediately before sowing. From each subplot, two samples were taken at a depth of 0-0.20 m and geo-referenced using the Leica Zeno 20 (Leica Geosystem, Heerbrugg, Canton St. Gallen, Switzerland). For each soil sample have been measured the content of sand (g kg −1 ), silt (g kg −1 ) and clay (g kg −1 ) by using the hydrometer method, pH with the pH meter method, Organic matter (g kg −1 ) with the Walkley-Black chromic acid wet oxidation method, the total nitrogen (g kg −1 ) with the Kjeldahl method and the C/N as the ratio of previous measures ( Table 2).

Experimental Design and Crop Management
The experimental site is a part of a long-term experiment (LTE) established in 1994 but still on-going [41] and consists of a rainfed 2-year rotation with durum wheat (Triticum turgidum L. subsp. durum cv. Grazia, ISEA) and maize (Zea mays L., DK440 hybrid Dekalb Monsanto, FAO Class 300) [17].
In this study, research activities focused on durum wheat plots where two soil management practices (main plot, 1500 m 2 ) and three nitrogen (N) levels (sub plot 500 m 2 ) were repeated in the same plots every year and arranged according to a split plot experimental design with two replications.
Plots with conventional tillage (CT), which is representative of the business as usual tillage practice in the study area, were ploughed along the maximum slope every year by a mouldboard plough (with 2 ploughs) at a depth of 0.40 m in autumn.
The seedbed was prepared with harrowing before the sowing date. The no tillage (NT) soil was left undisturbed and sprayed with herbicides before sowing prior to direct seed drilling.
The three N fertilizer treatments (N0, N1 and N2) corresponded to 0.90 and 180 kg N ha −1 were distributed in two rates. The unfertilized N0 treatment was chosen as control. The N1 treatment was compliant with the agri-environmental measures adopted within the Rural Development Plans at local scale. The N2 treatment was the business as usual N rates in the study area at the start of the experiment. Dates (dd/mm /yyyy) of all agronomic practices are reported in Table 3. Table 3. Agronomic management practices adopted during the two-year experimental period (dd/mm/yyyy).
At crop maturity (ZS92), the yield components have been measured by the number of kernels per spike (KS), the spikes m −2 and the grain yield (t ha −1 ) for both years.

Soil Plant Analysis Development (SPAD) Readings
For each plot and for each field survey day, three test areas were randomly selected and georeferenced the biomass sample positions by using the Leica Zeno 20 ( Figure 1) for a total of 36 ground control points (GCPs). At each test area, 10 fully expanded and intact leaves were chosen to make the readings with the SPAD (mod. Minolta 502). The readings were taken on the central and distal portion of the leaves, at the same time slot, around midday (11:00-13:00 a.m.).
The 10 leaves for each test area, after the SPAD readings, were cut and sealed in a plastic bag, placed in a portable refrigerator and transferred to the laboratory where the leaves' chlorophyll concentration (mg g −1 ) was measured.

Leaves Chlorophyll Concentration
The leaf chlorophyll concentration (mg g −1 ) was determined according to the Arnon's method [44]. The 10 fully expanded and intact acquired from the previous sampling activities and were weighed until we obtained about 0.1 g with a precision balance of α = 0.0001 g.
The leaf pigments were extracted in 80% acetone and leaf was pulverized completely by using the homogenizer Bio-Gen PRO200 (PRO Scientific Inc., Oxford, MS, USA) to obtain a completely white foliar tissue. Before the analysis, the leaf pigments were centrifuged for 5 min at 7 rpm to separate solid and liquid substance, and then more solution of 80% acetone was added until the achievement of a final volume of 25 mL. The determination of leaves chlorophyll concentration (mg g −1 ) is carried out by using a spectrophotometer Varian Cary 50 Scan ultraviolet (UV)-visible spectrophotometric (Agilent Technologies, Inc., Santa Clara, CA, USA), where the absorbance's (A) were measured at 663 and 645 nm. Total chlorophyll contents of each sample were computed from the follow equation Equation (1): mg total chlorophyll/g foliar tissue = 20.2 · (A645) + 8.02 · (A663) · V1000 · W (1) where: A = absorbance at specific wavelengths V = final volume of chlorophyll extract W = fresh weigh of tissue extracted.

Nitrogen Crop Status Determination
At each test area, fresh plant biomass was manually cut and collected in a 0.5 m long-row. The fresh plant biomass was oven-dried at 80 • C for 48 h and weighed (g). Before analyzing for the N content (%), the dry biomass was grounded to pass a 0.5 mm sieve.
The N content (%) was determined by the automated combustion analysis Dumas method [45,46] in an oxygen-enriched atmosphere at a high temperature by using the EA 1110 LECO CHNS-0 analyser (Leico Corporation, St. Joseph, MI, USA).
Starting from the N content (%) results, the NNI was calculated by dividing actual nitrogen concentration by the critical nitrogen concentration using the wheat dilution curve [47,48]. The following equations were used to calculate the NNI Equations (2) and (3): where: %N = actual N concentration N c = critical N concentration a = 5.35 DM = dry matter (g) b = 0.442.

Unmanned Aerial Vehicle (UAV) Flights, Image Acquisition and Processing
Multispectral images were acquired in correspondence of the three phenological stages (ZS22, ZS35 and ZS60).
Before each flight, a radiometric calibration with the calibration reflectance panel was performed. In addition, in order to perform better imagery acquisition, the UAV platform was equipped with the incident light sensor (ILS) which can record the variation of the light reflectance during the flight for each 9 spectral bands (wavelengths between 443 and 865 nm) of the MAIA Sentinel 2.
To obtain a high spatial accuracy of the remote-sensing images, 6 GCPs were georeferenced before each flight operation were used. All the relevant information regarding the flights is given in the following Table 4. Each image acquired by UAV flight operations required an image process to analyze all the relevant information. The image processing is composed by four main steps: raw to tiff conversion; orthomosaic reflectance map generation; selection and calculation of VIs; data extraction.
After exporting the raw images from the camera, the dedicated software MultiCam Stitcher Pro (SAL engineering and Eoptis, Russi, Italy) was used to convert the raw format image to the multi-channel tiff format.
To generate the orthomosaic reflectance maps the commercial software Agisoft Metashape (Agisoft LLC, St. Petersburg, Russia) was used. The software is based on structure from motion (SfM) algorithm [49]. The newly generated orthomosaic reflectance map was imported in R statistical software [50] and used to complete the remaining two main steps of the image processing.
To select the most relevant VI based on multispectral imagery for precision agriculture application, ten VIs were compared [51]. The VIs analysed in this study are reported in Table 5. Table 5. Vegetation index categories, calculated vegetation indices and related formula.

Vegetation Index
Formula References The csv file format contained the georeferenced positions of the biomass sampling was imported in R to extract the vegetation index values. To perform the values extraction, the "extract" function of "raster" R package [58] was used.

Yield Parameters
To characterize the yield for each of compared treatments, the number of kernels per spike (KS), the number of spikes m −2 and the grain yield (t ha −1 ) were measured at crop maturity (ZS92).
The KS was estimated on 30 spikes randomly collected per plot while the spikes m −2 was determined by counting the number of spikes along 1-m long row. The grain yield (t ha −1 ) expressed in dry matter, was measured by using a laboratory thresher for the three test areas (1 m long-row) per subplot.

Statistical Analysis
All statistical analyses were performed with R statistical software [50]. Every good analysis is based on a good model that is biologically relevant and based on realistic assumptions [59].
To select the model that can fit better the crop parameters and vegetation index dataset, three different models by using the "stats" [60] and "nlme" [61] R packages were built: (1) A one factor linear model was built by using the "lm" function on which the soil management was considered the main factor. (2) A full factorial linear model was also created by using the "lm" function on which the soil management, nitrogen input and year were considered as factors. (3) A mixed model was built by using "lme" function on which the soil management, nitrogen input and year were considered as fixed factors, and the Zadoks scale, the block and subplot were considered as random factor, according with [59].
In order to evaluate the best model, an 'a posteriori' selection was made, based on those statistics which put a penalty on 'complexity', such as the Akaike information criterion (AIC) [62] the Bayesian information criterion (BIC) [63], and likelihood ratio tests (LRTs) [64].
Based on the previous model analyses, the statistical model that fit better the crop parameters, vegetation index and the yield datasets was the mixed model, as reported by several studies [59,65].
Before performing the analysis of variance (ANOVA) it has been verified that the model met the three assumptions of the ANOVA. The Normality distribution of the model residual was checked both graphically (QQ-plot) and by performing the Shapiro-Wilk normality test. Moreover, the homoscedasticity was checked using the Levene test. The last ANOVA assumption was satisfied by the experimental design and the random sampling.
When all the three ANOVA assumptions were met, The ANOVA was applied to the model. Only when the ANOVA showed a significant difference between the factors (p-value < 0.05) the estimated marginal means post-hoc analysis was performed by using the "emmeans" function with the Bonferroni adjustment of the emmeans R package [66].
The linear relationships between the previous variables were analyzed, by using the Pearson correlation parameter (r) of the "cor" function of the R stats package [50], the coefficient of determination (R 2 ) and the root mean square error (RMSE).
To analyze the yield dataset the same statistical approach it was used as previously proposed for crop parameters and vegetation index dataset, but unlike before we set the year as a unique random factor while the soil management and nitrogen input were set as fixed factors.
The soil data were analyzed by using a full factorial linear model on which soil management and nitrogen input were considered as factors. The assumptions of the ANOVA application were checked and applied to the model. The investigation of the significant differences between groups factors were performed as before.

Crop Parameters
The ANOVA applied to the mixed model showed that the soil management (SM) was statistically significant for the NNI and SPAD readings while the nitrogen input (N) was statistically significant for each crop parameters (Table 6). *: Significant at p < 0.05%; **: Significant at p < 0.01%; ***: Significant at p < 0.001%.
The interaction between SM and N was statistically significant for all the measured variables (Table 6); this led to perform the estimated marginal means post hoc analysis on this factorial interaction.
The all pair-wise comparisons test showed that the factorial combination CT-N2 and NT-N2 performed statistically better than the other factorial combinations for the leaf chlorophyll concentration (mg g −1 ), %N, NNI and SPAD readings (Table 7) with a mean value of 2.73 (mg g −1 ), 2.18 (%), 1.19 and 47.94 respectively. The factorial combination CT-N1 and NT-N1 showed a statistical superiority over both unfertilized factorial combinations (CT-N0 and NT-N0) for each variable. This confirmed that confirmed that the nitrogen input is the key driver for the durum wheat growth development (Table 7). In detail, it was observed that there was an average percentage increase of +25.10%, +18.04%, +45.05%, +30.52% respectively for the leaf chlorophyll concentration (mg g −1 ), %N, NNI and SPAD readings.

Vegetation Index
The ANOVA applied to the mixed model showed that the SM was statistically significant only for 4 out of the 10 VIs analysed. These indices were the atmospherically resistant vegetation index (ARVI), green atmospherically resistant vegetation index (GARI), normalized difference red edge index (NDRE) and NDVI ( Table 8). The N factor had significant impact on all 10 VIs analysed with a p-value < 0.001. The interaction between SM and N was statistically significant only for the NIR-based VIs.
Those VIs analysed, which do not have the NIR band on the formula calculation, reported a statistical difference only between the N levels (Table 10). Both fertilized treatments showed a higher value than the unfertilized treatment for GLI and VDVI with a mean difference of +0.09 and +0.09 respectively. The VARI, VARI Green and Vig showed a superiority of N2, passing through the N1, to the N0 which showed the lowest values.

Regression Analysis
Both remote (through vegetation indices computation) and proximal sensing (through SPAD readings) data acquisition showed a significant linear relationship with leaves chlorophyll concentration (mg g −1 ) and NNI. The mean correlation parameter between all the remote and proximal sensing and leaves chlorophyll concentration (mg g −1 ) was R 2 = 0.47 and RMSE = 0.60. Furthermore, the average of the correlation parameters between the remote and proximal sensing with NNI was R 2 = 0.68 and RMSE = 0.31 (Table 11).
The MSAVI2 vegetation index showed a better correlation to be related to leaves chlorophyll concentration (mg g −1 ) with a R 2 value of 0.68 and RMSE equal to 0.47; for the NNI the NDRE vegetation index (Figure 2) showed a better performance with a R 2 of 0.85 and RMSE of 0.21. Table 11. Coefficient of determination and root mean square error (RMSE) of linear relationship between the leaves chlorophyll concentration (mg g −1 ), NNI and SPAD or vegetation indices.  The NIR-based VIs showed a better average correlation parameter for both crop parameters compared to the not NIR-based VIs (Table 12). The linear relationship between leaves chlorophyll concentration (mg g −1 ) and the VIs with the NIR band obtained a higher correlation parameter than the VIs without NIR band, with a difference of +0.32 for R 2 and −0.18 of RMSE (Table 12). A lower difference of +0.04 for R 2 and −0.03 RMSE for NNI was reported between VIs with NIR band and the VIs without the NIR band (Table 12).

Yield Parameters
The ANOVA showed a statistically significant impact of soil management (SM) on the number of kernels per spike (KS) and grain yield (t ha −1 ). The nitrogen input (N) affected significantly all the yield components and grain yield (t ha −1 ). The factorial combination SM × N showed a significant difference only on the KS and grain yield (t ha −1 ) (Table 13). With regard to the KS and grain yield (t ha −1 ) the estimated marginal means showed the same statistical difference as the VIs with the NIR band on the formula calculation, on which the NT-N2 and CT-N2 were statistically higher than the other factorial combinations (Table 14) with a mean value of 35 KS and 5.28 grain yield (t ha −1 ). The factorial combination CT-N1 and NT-N1 performed statistically better than both unfertilized factorial combinations, with a mean difference of +17 KS and +1.87 grain yield (t ha −1 ) ( Table 14).
The KS and grain yield (t ha −1 ) observed on the NT-N0 combination were statistically higher than the CT-N0 (Table 14), with a percentage increase of +41.18% and +41.83% respectively.
The number of spike m −2 reported a difference only between the three different nitrogen levels where the higher mean value was obtained by N2 with 567 spike m −2 , followed by N1 with 530 spike m −2 and N0 with 400 spike m −2 (Table 14).

Crop Parameters
Nitrogen input is the key driver for the durum wheat growth development [67]. As the nitrogen input increase, the overall crop nutritional status improves. This allowed a higher production of leaves' chlorophyll concentration (mg g −1 ), biomass nitrogen content (%), NNI and the SPAD readings as previously observed for durum wheat by [68] and [69].
On the unfertilized treatment, the NT performed better for each crop parameters than the CT. This is probably due to: the higher content of organic matter and consequently nitrogen availability present in the NT system 0-20 cm soil layer [70] than the compared CT plots ( Table 2).
The positive effect of CA compared with the conventional tillage practices was reported also by [71] who, by using the Decision Support System for Agrotechnology Transfer model (DSSAT) in the same experimental site, showed that NT leads to an increase in soil organic carbon with a rate of 0.43, 0.31 and 0.03 t ha −1 respectively for 180, 90 and 0 kg N ha −1 year −1 .
This characteristic of NT plots was observed also by [72] and [73] who recommended to adopt the CA to increase the soil organic matter, especially to mitigate the variations due to climate change [74].
Recently, [75] demonstrated that the use of the ARMOSA process-based crop model that CA can increase the annual soil organic carbon sequestration by 0.4% year −1 .

Vegetation Index
All VIs were susceptible to the different nitrogen input provided, therefore they are good indicators to identify crop variations due to different nitrogen input as observed by [76] where a significant difference was reported on the SPAD readings and NDVI values between 0, 80 and 160 kg N ha −1 . Moreover, [77] observed a significant difference on the NDVI values between 0, 55, 110 and 220 kg N ha −1 .
The ability of individual vegetation indices to detect different levels of nitrogen fertilization depends on the changes in reflectance recorded by each individual multispectral band. It is widely accepted that when a crop has a better overall nutritional status led to an higher leaves chlorophyll concentration (mg g −1 ) production [69]. A higher leaf chlorophyll concentration led to higher absorbance in the red and blue visible spectrum (VIS) area and reflect the green VIS area. This spectral characteristic of chlorophyll allowed the development of the Arnon method 1949 [44] with which we measured the leaves' chlorophyll concentration (mg g −1 ) and also the development of the optics that we used in this study such as the SPAD readings, multispectral camera and others proximal and remote-sensing tools applied in precision agriculture.
The VIs calculated with the NIR band detected the variation due to the SM within the N treatments, while the others VIs could not report a statistically alteration due the SM × N factorial interaction. The ARVI, MSAVI2, GARI, NDRE and NDVI statistically detected the superiority of NT-N0 over the CT-N0 as observed by [7] which in cotton observed a statistically increase of the NDVI on NT compared to CT.

Relationship Between SPAD Readings, Vegetation Indices and Crop Parameters
All significant linear relationships between the non-destructive methods and crop parameters were observed. It was reported by [78] that there was a significant linear relationship between 12 VIs and the nitrogen status in an experimental design where 8 different levels of N fertilization were applied on cotton (0, 130, 177, 194, 210, 307, 324 and 340 kg N ha −1 ). A significant and high linear relationship between 30 VIs and leaf chlorophyll content (LCC) was reported by [79], where the MDATT index reached a correlation accuracy of R 2 = 0.86 and a RMSE of 6.85.
The VI that reported the best correlation with the leaf chlorophyll concentration (mg g −1 ) was the MSAVI2, which was the same VI to obtain a significant relationship in each factorial combination reported by [68].
The NDRE was the VI to show the best relationship with NNI (Table 11), in accordance with [78], where for three different cotton growth stages NDRE was the best VI to be related with the nitrogen content.
The VIs which use the NIR band showed a higher mean correlation parameter with measured crop parameters than the VIs which do not use the NIR band on the formula calculation. This is probably due to the higher reflectance of the NIR in the plant cell, and in fact it is widely accepted that a plant with a good nutritional status has a higher reflectance in NIR spectral zone than a plant which does not have a good nutritional status [80].

Yield Components
The nitrogen input affected all the yield components, and this confirms what reported by [17] where the nitrogen was the key driver for the durum wheat growth and production after analyzing 20 years of data.
For the soil management the results found in the bibliography are discordant. [70] where only the unfertilized combination between CT and NT have been considered, reported the superiority of NT over the CT for KS and grain yield (t ha −1 ). While [17] reported that the NT did not provide any substantial advantage or disadvantage to durum wheat grain yield in comparison to CT. These results encourage in any case the adoption of the NT management because it leads to higher soil organic matter and nitrogen availability as described by [75] and allows farmers to have lower budget cost per ha as reported by [81].
The spike m −2 was affected by the different nitrogen levels as observed by [68]. The higher grain yield (t ha −1 ) obtained from NT-N0 over CT-N0 was +41.83% which is in line with the results obtained for KS (+41.18%).
The superiority of NT over CT is due to the higher organic matter (g kg −1 ). total nitrogen availability (g kg −1 ) and carbon nitrogen that NT plots had respect to CT ( Table 1), because of over 25 years of no soil-tillage practices [73].

Conclusions
Over 25 years of repeated no tillage practices allowed to NT plots to accumulate a higher soil organic matter content and total nitrogen availability than the CT plots.
In order to obtain valid results from a split plot LTE and interpret them correctly, the use of the mixed model is recommended, because it is the statistical model capable of better fitting the experimental data than the classics linear models.
Nitrogen input is a key driver for durum wheat growth development and production, as the nitrogen input increases. There is an increase also in the crop parameters analyzed. When no mineral nitrogen was supplied to the crop, NT performed better than CT for any parameters measured in this study apart from the spike m −2 .
All VIs were capable of detecting variation due to the nitrogen input, but few of them detected the difference due to the interaction between soil management and nitrogen input, which were the VIs with the NIR band on the formula calculation. Moreover, we recommend readers use the MSAVI2 to monitor the leaf chlorophyll concentration while we suggest using the NDRE to monitor the nitrogen durum wheat status.
The presented results can also provide an important contribution in the development and implementation of agro-environmental policies aimed at containing production inputs to reduce the environmental impact of cereal-based cropping systems.

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