Detection of Spatial and Temporal Variability of Wheat Cultivars by High-Resolution Vegetation Indices

An on-farm research study was carried out on two small-plots cultivated with two cultivars of durum wheat (Odisseo and Ariosto). The paper presents a theoretical approach for investigating frequency vegetation indices (VIs) in different areas of the experimental plot for early detection of agronomic spatial variability. Four flights were carried out with an unmanned aerial vehicle (UAV) to calculate high-resolution normalized difference vegetation index (NDVI) and optimized soil-adjusted vegetation index (OSAVI) images. Ground agronomic data (biomass, leaf area index (LAI), spikes, plant height, and yield) have been linked to the vegetation indices (VIs) at different growth stages. Regression coefficients of all samplings data were highly significant for both the cultivars and VIs at anthesis and tillering stage. At harvest, the whole plot (W) data were analyzed and compared with two sub-areas characterized by high agronomic performance (H) yield 20% higher than the whole plot, and low performances (L), about 20% lower of yield related to the whole plot). The whole plot and two sub-areas were analyzed backward in time comparing the VIs frequency curves. At anthesis, more than 75% of the surface of H sub-areas showed a VIs value higher than the L sub-plot. The differences were evident also at the tillering and seedling stages, when the 75% (third percentile) of VIs H data was over the 50% (second percentile) of the W curve and over the 25% (first percentile) of L sub-plot. The use of high-resolution images for analyzing the frequency value of VIs in different areas can be a useful approach for the detection of agronomic constraints for precision agriculture purposes.


Introduction
Monitoring the spatial and temporal variability of wheat within a season is crucial to decision-making in precision farming. Precision agriculture is a modern farming management concept using digital techniques to monitor and optimize agricultural production processes. Precision agriculture can play an important role in enhancing crop yield and ensuring sustainability [1]. Among the tools used to acquire information, unmanned aerial vehicles (UAVs) equipped with visible and near-infrared cameras, provide, in a fast and easy way, field data for precision agriculture applications [2,3]. The resolution of information from satellite data typically ranges from 5 to 30 m pixels and is unsuitable in agronomy trials given the limitations of real-time monitoring and accuracy [4]. In contrast to satellite imagery and aircraft-based remote sensing, UAVs can be used frequently during the entire growth period [5]. Furthermore, vegetation indices (VIs) of UAV imagery have the same ability as ground-based recordings to quantify crop responses to experimental treatments [6]. UAVs are a useful technology for crop monitoring at different scales and can be used for agronomic experiments was harvested in July 2015. The soil was prepared with the help of minimum-tillage equipment and basic fertilization was performed with 200 kg ha −1 of P 2 O 5 (27%). Nitrogen was applied at the end of March (tillering) with a slow release nitrogen fertilizer, 35% N; 23% SO 3 , and at April (booting/stem elongation) with 150 kg ha −1 of fertilizer with 30% N; 23% SO 3 and with two foliar nitrogen application. The weed and pest control scouting was carried out with chemical pesticides during the crop cycle.
A standard agro-meteorological station was placed in the experimental fields and temperatures, rainfall, wind, humidity, and radiation were recorded. The minimum temperature (−1 • C) was recorded in December 2014, while the maximum (36.2 • C) was recorded in July 2015. Through the whole cropping season (December 2014-July 2015), precipitation amounted to 680.2 mm; a good half of the total precipitation was recorded in 6 weeks, from the final third of January to the first third of March. From tillering to anthesis, rainfall was scarce, while from anthesis to ripening, rainfall was considerable. Air temperatures were mild from emergence to anthesis and afterwards, maximum air temperatures were around 30 • C, as expected.
The phenological stages of the wheat crop were periodically recorded according to the Zadoks Scale (ZS), which is a standardized reference scale used to evaluate and measure plant growth stage in cereals [25]. During the crop growth, sampling was carried out for each cultivar: 10 at seedling growth 13-ZS (February 23), 12 at tillering 25-ZS (March 30), 11 at anthesis 65-ZS (May 14), and 11 at harvest-99 ZS (July 7). Plants from 0.5 m 2 were georeferenced and hand cut for the calculation of yield-related traits (biomass, green leaf area index-LAI, number of spikes, plant height, yield, and thousand kernels weight). For each cultivar, 44 georeferenced samplings during crop growth were collected.
Whole plant dry mass was determined after oven drying the fresh plant material at 75 • C until constant weight for 48 h. The sampling points based on visual-spatial crop variability and identified by low-resolution orthoimages of flight were processed in real time in the field. Ground truth coordinates of target locations were recorded with GPS Leica Viva GS15 (Leica Geosystems AG, Heerbrugg, Switzerland) at each sampling and flight.
Furthermore, at harvest, within each cultivar, yield data and yield components of the whole plot (720 m 2 ) were measured, as well as for four sub-areas, 10 m 2 each was characterized by the different yield levels in Figure 1 (H = high yield and L = low yield). The sub-areas were selected within the homogeneous areas detected using a hierarchical clustering method, Ward's minimum variance approach, as reported in a previous paper for Odisseo [25] and calculated for Ariosto. NDVI and OSAVI high-resolution data of the four areas were analyzed backward in time to deliver relative frequency histograms of VI values. applied at the end of March (tillering) with a slow release nitrogen fertilizer, 35% N; 23% SO3, and at April (booting/stem elongation) with 150 kg ha −1 of fertilizer with 30% N; 23% SO3 and with two foliar nitrogen application. The weed and pest control scouting was carried out with chemical pesticides during the crop cycle. A standard agro-meteorological station was placed in the experimental fields and temperatures, rainfall, wind, humidity, and radiation were recorded. The minimum temperature (−1 °C) was recorded in December 2014, while the maximum (36.2 °C) was recorded in July 2015. Through the whole cropping season (December 2014-July 2015), precipitation amounted to 680.2 mm; a good half of the total precipitation was recorded in 6 weeks, from the final third of January to the first third of March. From tillering to anthesis, rainfall was scarce, while from anthesis to ripening, rainfall was considerable. Air temperatures were mild from emergence to anthesis and afterwards, maximum air temperatures were around 30 °C, as expected.
The phenological stages of the wheat crop were periodically recorded according to the Zadoks Scale (ZS), which is a standardized reference scale used to evaluate and measure plant growth stage in cereals [25]. During the crop growth, sampling was carried out for each cultivar: 10 at seedling growth 13-ZS (February 23), 12 at tillering 25-ZS (March 30), 11 at anthesis 65-ZS (May 14), and 11 at harvest-99 ZS (July 7). Plants from 0.5 m 2 were georeferenced and hand cut for the calculation of yield-related traits (biomass, green leaf area index-LAI, number of spikes, plant height, yield, and thousand kernels weight). For each cultivar, 44 georeferenced samplings during crop growth were collected.
Whole plant dry mass was determined after oven drying the fresh plant material at 75 °C until constant weight for 48 h. The sampling points based on visual-spatial crop variability and identified by low-resolution orthoimages of flight were processed in real time in the field. Ground truth coordinates of target locations were recorded with GPS Leica Viva GS15 (Leica Geosystems AG, Heerbrugg, Switzerland) at each sampling and flight.
Furthermore, at harvest, within each cultivar, yield data and yield components of the whole plot (720 m 2 ) were measured, as well as for four sub-areas, 10 m 2 each was characterized by the different yield levels in Figure 1 (H = high yield and L = low yield). The sub-areas were selected within the homogeneous areas detected using a hierarchical clustering method, Ward's minimum variance approach, as reported in a previous paper for Odisseo [25] and calculated for Ariosto. NDVI and OSAVI high-resolution data of the four areas were analyzed backward in time to deliver relative frequency histograms of VI values.

UAV System and Flight Missions
Four flight missions were carried out in the same dates to measure the destructive plant samples: for seedling (February 23), tillering (March 30), anthesis (May 14), harvest (July 7), with eBee UAV (senseFly, Cheseaux-sur-Lausanne, Switzerland). The eBee, designed as a fixed-wing UAV for application in precision agriculture (payload of 150 g). All flights were carried out in stable ambient light conditions from 11:30 a.m. to 1:00 p.m., with excellent visibility and wind below 5 m s −1 , at a flight altitude of 97 m (A.G.L.). The imaged area of the experimental field, including the surroundings, is about 15 hectares so that 96 camera stations (single flight) and 15 min flight were needed. The autopilot analyzes (continuously) the inertial measurement unit (IMU) and onboard GPS data to control every aspect of the eBee's flight. The integration of the UAV and sensors with GPS and IMU enables obtaining direct georeferencing imaging after image processing. At any acquisition date, two flights were carried out, the first one by a Canon Powershot S110 photo camera (visible spectrum, RGB-red/green/blue) for visible RGB image (orthophoto) to run a rapid analysis for visual crop variability. The second flight used a Canon Powershot S110 NIR camera (near infra-red, NIRGB-near infra-red/green/blue) which provides maximum absorption peaks at 550 nm (green), 625 nm (red), and 850 nm (NIR) wavelengths, respectively, allowing the computation of VIs.
To avoid geometric distortion due to low altitude, 96 overlapping pictures from each camera and flight were used for mosaicking to produce an ortho-image. An 80% frontal overlap and an 80% side overlap were used.
To relate and orient UAV imagery to the ground, ten ground control points (GCPs) were distributed across the field at the beginning of the season, to obtain photogrammetric imagery with uniform vertical and horizontal accuracy. The GCPs were 25 cm × 25 cm square, with a specific albedo for camera calibration (atmospheric corrections). The GCP coordinates were ensured with a Leica Viva GS15 (Leica Geosystems AG, Heerbrugg, Switzerland) GPS (horizontal accuracy of 0.025 m-vertical accuracy of 0.035 m). Fixed targets were used for a more accurate geo-referencing of UAV aerial imagery and for overlaying the measurements from multiple dates.

Data Processing
The acquired images were processed by the eMotion software (senseFly, Cheseaux-sur-Lausanne, Switzerland) to generate low-resolution visible (RGB) images of the crops. The eBee's supplied software to build a project using the drone's geo-tagged images. The project is loaded on a laptop in Pix4Dmapper Ag (senseFly, Cheseaux-sur-Lausanne, Switzerland) to run a rapid analysis for visual crop variability.
The multiple overlapped images of the whole plot were stitched and ortho-rectified to create the geo-referenced ortho-mosaicked image. In the laboratory, data processing (orthomosaicking) of acquired images were performed with Postflight Terra 3D software package, a customized version of the Pix4D digital photogrammetric solution specifically optimized for the eBee, to generate ortho-images. Postflight Terra 3D incorporates a scale-invariant feature transform (SIFT) algorithm to match key points across multiple images [26] and process data in three key steps: (1) initial processing, (2) point cloud densification, and (3) orthomosaic generation.
Ortho-rectification by aero-triangulation and mosaicking were elaborated for processing, starting from the exterior position and orientation parameters provided by the UAV internal system (roll, pitch, and yaw angles). Orthoimages were produced from the flights, with a pixel resolution of 5 cm. Ten GCPs (Leica Viva GS15 (Leica Geosystems AG, Heerbrugg, Switzerland) were used to complete the external camera orientation. The orthomosaic was georeferenced to UTM-WGS84 zone 33N Italy. The final outputs were an RGB (Visible) GeoTIFF with a resolution of 5 cm 2 pixels. The normalized difference vegetation index (NDVI) and the soil-adjusted vegetation index (SAVI) layers were generated in a raster calculator from extracted red (R) and near infra-red (NIR) channels. The index calculator function of Postflight Terra 3D was used for generating VI maps. A certificate of calibration of Canon S110 NIR was uploaded in the software to optimize internal camera parameters, such as focal length, principal points, and lens distortions. The ten GCPs with a known albedo for Red, Green, and NIR channel (reflectance panel) were used to calibrate the camera to achieve uniform quality of image (exposure and brightness) and for atmospheric correction in the software section "processing options," point 3 DSM, orthomosaic, index and for creating VIs map. The resolution of the reflectance map has been set at 5 cm 2 pixels, the whole map was elaborated for each VI and date, and the NDVI and OSAVI formula were selected (index map function). The final outputs were an NDVI GeoTIFF and an OSAVI GeoTIFF, both with a resolution of 5 cm 2 pixels. GeoTIFF images and georeferenced sampling data were processed for agronomic purpose with QGIS 2.8.1.

Vegetation Indices
The normalized difference vegetation index (NDVI) was calculated according to Equation (1): The NDVI ranges from −1.0 to 1.0, where positive values indicate increasing greenness and negative values indicate non-vegetated features. It has some disadvantages though such as saturation in later growth stages [27,28].
The optimized soil-adjusted vegetation index (OSAVI) was calculated according to Equation (2). The OSAVI index was proposed by Rondeaux et al. [29] using reflectance in the red and near-infrared bands; through the following formulation: where 0.16 is the soil adjustment coefficient, selected according to the optimal value to minimize soil background variations.

Statistical Analysis
Regression analysis, coefficients of determination, significance levels, and RMSE were computed on two sets of geo-referred data: ground crop samplings (LAI, Biomass, and spikes for square meter) and remote UAV data using the statistical package Origin PRO 8 (Origin Lab Corporation, Northampton, MA, USA). H and L sub-area were selected starting from homogeneous areas identified in a previous paper [22] as follow: the yield-related traits and VIs data were analyzed using the hierarchical clustering Ward's minimum variance approach [30] to classify observations into groups, in which the group members have common properties. Statistical procedures were computed using OriginPRO 8. The difference between H and L sub-areas were recorded by the analysis of variance.

Results and Discussion
Ground agronomic data (biomass, green leaf area index, spikes, and plant height) of each cultivar have been correlated to high-resolution multispectral images (NDVI and OSAVI). The correlation of the whole crop cycle data was performed in accordance with van Ittersum et al. [31], who states that it is essential to study the relationship among VIs and crop traits, to evaluate and estimate the yield potential and the yield gap (Table 1). In the present study, polynomial regressions were found to be the simplest adjustment to report the time-course variability along the growth stage, starting from seedling to harvest. Moreover, significant relationships with different R 2 values were found for both the indices and crop growth traits. The spectro-radiometric response of the wheat canopy was not linear during the growing period in accordance with Aparicio et al. [32] and Dang et al. [33], since the high sensitivity of the VIs when LAI is lower (early season) and less sensitive after the canopy closes.
The relationships among VIs and crop parameters were highly significant: the R 2 for NDVI vs. biomass reached 0.953 for Odisseo and 0.859 for Ariosto, while the R 2 of OSAVI vs. biomass was 0.943 for Odisseo and 0.857 for Ariosto. The relationships among each VI and LAI showed an R 2 of 0.910 for NDVI Odisseo and R 2 of 0.777 for Ariosto, and it was 0.879 for OSAVI Odisseo and 0.763 for Ariosto. The relationships among each VI and height showed an R 2 of 0.828 for NDVI Odisseo and 0.817 for OSAVI Odisseo and R 2 of 0.518 for both NDVI and OSAVI for Ariosto. Ground truth georeferenced data of the number of spikes per square meter at harvest stages were also related to NDVI and OSAVI data; significant relationships were found with an R 2 of 0.788 for Odisseo and 0.769 for Ariosto recorded for NDVI index and an R 2 of 0.790 for Odisseo and of 0.597 for Ariosto recorded for OSAVI index.
Briefly, all R 2 of yield components were highly correlated (p-value < 0.001) with VIs, and are higher for Odisseo than the Ariosto.

Agronomic Data and VIs at Different Crop Growth Stages
In the whole plot, the biomass at seedling stage ranged from 68 to 156 g m −2 for both cultivars and that of LAI (0.12 to 0.47) ( Table 2), while VIs ranged from 0.26 to 0.57 for NDVI and from 0.38 to 0.65 for OSAVI, for both cultivars. The light differences recorded within yield components as well as VIs value gave no significant relationships among VIs and yield components (Table 3). This was due to a low amount of accumulated biomass and LAI values, as well as the difficulties of simple indices to discriminate between soil and vegetation as found also by many authors, such as Aparicio et al. [34] and Chlingaryan et al. [24].
At the tillering stage, the spatial variability of crop traits was highly evident; thus, the most productive data was shown by the biomass and LAI which was four times higher than the less favored zones ( Table 2).
In the most productive zones, VI values ranged from 0.85 to 0.99 (for both cultivars), which revealed 2.5 times more than the lower value ( Figure 2).
The relationships between VIs and crop parameters were highly significant for both cultivars (Table 3). For biomass, the R 2 reached 0.88 for Odisseo and 0.66 for Ariosto for LAI, the R 2 reached 0.85 for Odisseo and 0.69 for Ariosto. The plant height and VI regression were not significant for Ariosto. These results are in accordance with findings of several authors at tillering, for example Dalla Marta et al. [35] found a highly significant (positive) relationship with biomass weight, while Reyniers and Vrindts [36] and Magney et al. [37] found respectively an R 2 of 0.76 and R 2 = 0.62 between NDVI and biomass.
At anthesis, the high spatial variability of the crop traits value and VI values were still evident for both cultivars (Table 2, Figure 2). VIs were highly correlated to crop traits as reported in Table 3; the values of the R 2 for NDVI vs. biomass reached 0.89 for Odisseo and 0.82 for Ariosto, while the R 2 of OSAVI vs. biomass was 0.80 for Odisseo and 0.78 for Ariosto (p < 0.001).
Villegas et al. [39] that reported the highest values of R 2 for the relationship between NDVI and crop dry weight at anthesis. Cabrera-Bosquet et al. [40] have demonstrated a strong linear regression between durum wheat and VIs. On the contrary, Dalla Marta et al. [35] found completely different results, with no correlations observed between the crop parameters and the indices, due to VIs saturation.  Table 2. Range value (min and max), mean and standard deviation (SD) of yield components (biomass, LAI, plant height, and spike) and vegetation indices (VIs) (NDVI and OSAVI) measured at seedling, tillering, and anthesis stages.    These results are in accordance with the Marti et al. [38], who found the highest correlation values between growth and NDVI measurements when performed around anthesis, such as with Villegas et al. [39] that reported the highest values of R 2 for the relationship between NDVI and crop dry weight at anthesis. Cabrera-Bosquet et al. [40] have demonstrated a strong linear regression between durum wheat and VIs. On the contrary, Dalla Marta et al. [35] found completely different results, with no correlations observed between the crop parameters and the indices, due to VIs saturation. Table 4 reports yield and yield components at harvest for the complete plot (W), as well as data of the two sub-areas, H O and H A with high yield and L O and L A with low yield. Table 4. Yield (t ha −1 ), biomass (g m −2 ), spike number (n • m −2 ), and plant height (cm) of Odisseo and Ariosto at harvest stage. Mean data are related to the whole plot (W) and to the sub-areas; H o (high yield of Odisseo) and H A (high yield of Ariosto), L O (low yield Odisseo), and L A (low yield Ariosto). The p-level values are related to H and L significant differences within each cultivar. The yield data of both cultivars were in accordance with those reported by Visioli et al.

W H O (SD) L O (SD) p-level W H A (SD) L A (SD) p-level
[41] in Italy. The crop yield of the whole plots was 5.73 t ha −1 for Odisseo (W O ) and 5.05 t ha −1 for Ariosto (W A ). The yield of H O and H A were about 20% higher than the whole plot yield, which in turn were 20% higher than L O and L A . Yield components were ranked according to the yield levels: biomass and spikes per square meter of H were about 25-35% higher than the W and 37-50% higher than L sub-plot areas . Differences in plant heights were less marked. The subplot detected at harvest was analyzed backward in time from anthesis to seedling. At anthesis, the analysis of VI value frequencies showed the maximum evidence of different crop status of H and L sub-areas and W plots. More than the 75% of the surface of the sub-areas showed a value higher (H sub-plot) or lower (L sub-plot) than W plot (Figure 3). The H and L areas curves showed a low overlapping. The H curve appeared narrow and tall. Data were confirmed by all the yield components that were significantly different from each other in the two sub-areas and W plot.
At tillering, the histogram of VIs frequency (Figure 4) shows differences among the whole plot data, both H and L areas. These differences were lower than those detected at the anthesis stage.
For both the cultivars and indices, the H zones showed values of the frequency which were higher than the 75% of H O and H A , and L A , L O VI values were higher and lower respectively than the median value of the whole plot. In all graphs, about 75% of the H data (third percentile) were higher than the second percentile of the W curve as well as the first percentile of L data. The highest differences between H and L curves were detected for the OSAVI index, Ariosto (d), because the index that considers soil brightness are perform slightly better than the NDVI, as reported in different studies [42,43]. Table 4. Yield (t ha −1 ), biomass (g m −2 ), spike number (n° m −2 ), and plant height (cm) of Odisseo and Ariosto at harvest stage. Mean data are related to the whole plot (W) and to the sub-areas; Ho (high yield of Odisseo) and HA (high yield of Ariosto), LO (low yield Odisseo), and LA (low yield Ariosto). The p-level values are related to H and L significant differences within each cultivar.    In all graphs, about 75% of the H data (third percentile) were higher than the second percentile of the W curve as well as the first percentile of L data. The highest differences between H and L curves were detected for the OSAVI index, Ariosto (d), because the index that considers soil brightness are perform slightly better than the NDVI, as reported in different studies [42,43].    At the seedling stage, the relative frequency curves of the two VIs are presented in Figure 5. In all the graphs, the curves appeared to be very close, except for the OSAVI of Ariosto (d) at the seedling stage.

Odisseo
In all graphs, about 75% of the H data (third percentile) were higher than the second percentile of the W curve as well as the first percentile of L data. The highest differences between H and L curves were detected for the OSAVI index, Ariosto (d), because the index that considers soil brightness are perform slightly better than the NDVI, as reported in different studies [42,43].

Conclusions
In the last few years, UAVs have been used as tools that provide a high number of information per square meter for precision agriculture management and its development. The potentials of the information provided by the drone images have not yet been fully explored. In this paper, we have analyzed the frequency distribution curves of high-resolution VIs images of the two cultivars of durum wheat.
Within each cultivar, two yield sub-areas were detected by cluster analysis at the harvest stage and analyzed backward at different growth stages, in comparison with the whole plot data. The crop yield of the whole plots revealed 5.73 t ha −1 for Odisseo (WO) and 5.05 t ha −1 for Ariosto (WA). The yield for HO and HA recorded 20% higher than the whole plot yield, which in turn were 20% higher than LO and LA, moreover the biomass and the spikes per square meter of H was noted 25-35% higher than the W and 37-50% higher than L sub-plot areas, confirming the ability of cluster analysis to detect groups with the same variability.
The analysis of frequencies of high-resolution VI images showed that the 75% of H data (third percentile) were recorded higher than the second percentile of the W curve as well as for the first percentile of L data. The maximum evidence of different crop status of H and L sub-areas and W plots were detected at anthesis by both indices and cultivars.
The analysis of frequency curves provided significant differences at the seedling stage where agronomic parameters showed only slight differences which were revealed as not significant following the traditional approach (regression between agronomic traits and VIs). For Ariosto, the seedling OSAVI index showed an interesting differentiation among areas, confirming the ability of OSAVI to consider the soil brightness.
The analysis of frequencies of high-resolution VI images provided detailed information on crop spatial and temporal variability and can allow developing a new approach for automatic detection of suitable vegetation indices, and spatial and temporal crop variability for precision agriculture practices.
Author Contributions: The authors, S.M. and A.A. contributed equally to this work.

Conclusions
In the last few years, UAVs have been used as tools that provide a high number of information per square meter for precision agriculture management and its development. The potentials of the information provided by the drone images have not yet been fully explored. In this paper, we have analyzed the frequency distribution curves of high-resolution VIs images of the two cultivars of durum wheat.
Within each cultivar, two yield sub-areas were detected by cluster analysis at the harvest stage and analyzed backward at different growth stages, in comparison with the whole plot data. The crop yield of the whole plots revealed 5.73 t ha −1 for Odisseo (W O ) and 5.05 t ha −1 for Ariosto (W A ). The yield for H O and H A recorded 20% higher than the whole plot yield, which in turn were 20% higher than L O and L A , moreover the biomass and the spikes per square meter of H was noted 25-35% higher than the W and 37-50% higher than L sub-plot areas, confirming the ability of cluster analysis to detect groups with the same variability.
The analysis of frequencies of high-resolution VI images showed that the 75% of H data (third percentile) were recorded higher than the second percentile of the W curve as well as for the first percentile of L data. The maximum evidence of different crop status of H and L sub-areas and W plots were detected at anthesis by both indices and cultivars.
The analysis of frequency curves provided significant differences at the seedling stage where agronomic parameters showed only slight differences which were revealed as not significant following the traditional approach (regression between agronomic traits and VIs). For Ariosto, the seedling OSAVI index showed an interesting differentiation among areas, confirming the ability of OSAVI to consider the soil brightness.
The analysis of frequencies of high-resolution VI images provided detailed information on crop spatial and temporal variability and can allow developing a new approach for automatic detection of suitable vegetation indices, and spatial and temporal crop variability for precision agriculture practices.
Author Contributions: The authors, S.M. and A.A. contributed equally to this work.