1. 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.
2. Materials and Methods
2.1. 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.
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).
2.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 m2) and three nitrogen (N) levels (sub plot 500 m2) 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.
2.3. Measurements
In both growing seasons under analysis (2018 and 2019) have been acquired Soil Plant Analysis Development (SPAD) readings (Konica Minolta Sensing 2003, Osaka, Japan), leaves chlorophyll concentration (mg g
−1), N content (%), nitrogen nutrition index (NNI) at the crop tillering, stem elongation and anthesis phenological stages (ZS22, ZS35 and ZS60 respectively were ZS = Zadoks scale [
43]).
In addition, multispectral images have been acquired (MAIA S-2, SAL engineering, Russi, Italy) by using an unmanned aerial vehicle (UAV) platform (DJI Matrice 600 pro, SZ DJI Technology Co., Ltd., Shenzen, China) to compute several vegetation index (VI) algorithms.
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.
2.3.1. 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.
2.3.2. 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):
where:
A = absorbance at specific wavelengths
V = final volume of chlorophyll extract
W = fresh weigh of tissue extracted.
2.3.3. 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:
2.3.4. 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 geo-referenced 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.
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.
2.3.5. 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.
2.4. 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.
5. 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.