An Integrated Approach to Assessing the Soil Quality and Nutritional Status of Large and Long-Term Cultivated Rice Agro-Ecosystems

The aim of this study is to develop an integrated approach to soil quality and fertility assessment in high-yielding rice agro-ecosystems threatened due to overexploitation of soil resources by intensive agriculture. The proposed approach is implemented considering representative pilot fields allocated throughout a study area based on the assumption that soils of similar general properties present a similar nutritional status due to common long-term management practices. The analysis includes (a) object-based image analysis for land zonation, (b) hot-spot analysis for sampling scheme evaluation, (c) setting of critical thresholds in soil parameters for detecting nutrient deficiencies and soil quality problems, and (d) Redundancy Analysis, TITAN analysis, and multiple regression for identifying individual or combined effects of general soil properties (e.g., organic matter, soil texture, pH, salinity) or non-soil parameters (e.g., topographic parameters) on soil nutrients. The approach was applied using as a case study the large rice agro-ecosystem of Thessaloniki plain in Greece considering some site specificities (e.g., high rice yields, calcareous soils) when setting the critical thresholds in soil parameters. The results showed that (a) 62.5% of the pilot fields’ coverage has a simultaneous deficiency in Zn, Mn, and B, (b) organic matter (OM) was the most significant descriptor of nutrients’ variance, and its cold spots (clustered regions of low OM values) showed important overlapping with the cold spots of K, Mg, Zn, Mn, Cu, and B, (c) a higher rate of availability increase in P, K, Mg, Mn, Zn, Fe, Cu, and B was observed when the OM ranged between 2 and 3%, and (d) the multiple regression models that assess K and P concentrations based on general soil properties showed an adequate performance, allowing their use for general assessment of their soil concentrations in the fields of the whole agro-ecosystem.


Introduction
The green revolution and the continuous introduction of advanced agricultural equipment and genetically improved rice varieties for achieving higher yields has substantially increased the pressures on natural resources, such as water and soil [1][2][3].Despite the fact that the benefits of these actions are especially highlighted through cereals, which have shown more than double mean global yields since the 1960s [4], there are concerns about the viability or the prospect of further yield increase, not only due to uncertainties in water availability and climate change [5,6], but also due to soil quality degradation [7] or combined nutrient and water limitations [8].
The high probability of soil quality degradation due to intensive agriculture may inhibit the increase of yields, and this had been already indicated more than two decades ago by Cassman and Pingali [9,10], which presented long-term experiments (20-25 years) of irrigated rice systems in different environments and soils that started during the 1960s and 1970s.Even though the rice varieties used in their experiments were regularly replaced with the highest-yielding and the most resistant germplasm to diseases and insects available at each point in time, negative yield trends were observed, probably due to soil quality degradation.Despite the fact that many years have passed, the title of their paper "Intensification of Irrigated Rice Systems: Learning from the Past to Meet Future Challenges" seems to be more than up-to-date not only for rice but also for all crops, since the race for increasing yields has set the target of doubling global crop production by 2050 in order to meet the projected demands from the rising population [11].Unfortunately, the current yield trends are insufficient to achieve this goal [11], and there is general evidence that the relative yield gain decreases over time, while in many cases yield plateaus or abrupt decreases in the rate of yield gain have been observed, including rice in eastern Asia and wheat in Northwest Europe, which share a very large portion of total global rice and wheat production [12].
The case of irrigated rice probably constitutes the best paradigm to highlight the aforementioned issues, because: • the flooding conditions have a strong impact on physical [13][14][15] and chemical soil properties [16] with evident morphological/pedological alterations after long-term rice cultivation [17] in comparison to non-flooded adjacent territories with other crops; and • rice was the first crop to be completely sequenced [18], providing availability of high-quality genomes and long evolutionary trajectories that contain a virtually untapped reservoir of genes and traits that can be used for breeding a new generation of sustainable rice cultivars, known as Green Super Rice [19].The new era of rice varieties may reduce the risk of yield failure due to adverse environmental conditions but at the same time could be the reason for soil degradation through excessive nutrient depletion by high yields.
Apart from the evidence provided by Cassman and Pingali [9,10] about rice crop pressures on soil quality and productivity, Tong et al. [20] also provided more striking evidence, showing that the rice yield production per unit chemical fertilizer decreased from 164 kg kg −1 in 1961 to 10 kg kg −1 in 1998 in China.This probably occurred due to the fact that N use rates for 45% of rice paddy fields in China are 20-40% higher than the environmentally optimal N rates [21].Ali [22] analyzed land degradation due to technological change and cultivation of high-yielding variety (HYV) rice, using power tillers, low-lift irrigation pumps, chemical fertilizers, and pesticides during the period 1985-2000 in southwestern Bangladesh.The analysis showed significant changes in soil texture, bulk density, and structure, a soil salinity increase, and depletion of macronutrients and organic C in double-cropped irrigated rice fields, while power-tilled plots showed a higher rate of soil erosion and loss of soil structure than the traditionally ploughed fields.Soil degradation affected the yield of rice, while cultivation of HYV rice decreased the fodder production, affecting the livestock population.A similar type of soil quality degradation was also observed in rice-shrimp cropping systems of the same territory [23].Lima et al. [24] observed that land degradation and, especially, physical soil degradation of rice soils in South Brazil is the main obstacle for growing other crops in rotation with rice.
The available data on rice from Food and Agriculture Organization of the United Nations (FAO) [4] of the period 1961-2016 for the three most productive countries in terms of yield in the Mediterranean area (Egypt, Turkey, and Greece) highlight the yield increase that succeeded after 1960 (Figure 1).These three countries are also in the group of the five most productive countries in terms of yield in the world, together with Australia and the United States.Egypt was the first country in the world that passed the threshold of 10 Mg ha −1 at the country level during 2006 (Figure 1); however, since then, the yields showed a declining trend during the period 2006-2016, which created a lot of questions about the possible reasons for the decline.According to a commodity intelligence report of United intelligence report of United States Department of Agriculture (USDA) [25], the yields in 2015 in Egypt were lower than in the previous year due to extreme heat, but further justifications for the general declining trend of the period 2006-2016 have not been documented.A similar drop in yield was also observed in the case of Turkey after passing the threshold of 9 Mg ha −1 at the country level during 2011, while Greece presents a concave yield plateau during 1995-2016 (Figure 1).The aforementioned cases of Mediterranean countries show that a further rice yield increase at the country level will be difficult to obtain, and special attention should be given not only to the consequences of climate change but also to possible impacts on soil resources.The primary and most important procedure for assessing soil quality status is soil analysis.Zonation through remote sensing techniques (RSTs) and geographical information systems (GISs) have played a key role for a) optimizing the number of sampling sites for the development of representative sampling networks based on the detection of spatial soil differences and crop growth anomalies [26][27][28], b) the identification of the possible origin of such differences/anomalies through in situ observations/samplings, and c) the design of faster and more targeted interventions.On the other hand, it is very important to stress the high cost of such actions in large agricultural areas and the need for additional indirect methods in order to further reduce their cost through even more targeted soil samplings and soil management interventions.
The aim of this study is to develop an integrated approach to soil quality and fertility assessment in high-yielding rice agro-ecosystems threatened due to overexploitation of soil resources by intensive agriculture.Considering the urgent need for soil monitoring, quality evaluation, and management in rice agro-ecosystems, this approach could contribute to decision-making by farmers and management authorities.
A variety of methodologies, including remote sensing, geostatistics, and multiregression analysis were applied as nodes of an integrated approach.The overall approach was applied to extended pilot fields assuming that soils of similar general properties present a similar nutritional status in large agro-ecosystems where common cropping management practices are followed for a long-term period.The robustness of the aforementioned assumption is indirectly evaluated through measures of statistical significance that derive from the proposed methodologies.The proposed approach is presented in this study using as a case study the large rice agro-ecosystem of Thessaloniki plain in Greece.

Study Area
The study area is the agricultural plain downstream of Axios river in Greece (Figure 2).The mean annual precipitation and temperature in the study area are 440 mm and 15.8 o C, respectively, The primary and most important procedure for assessing soil quality status is soil analysis.Zonation through remote sensing techniques (RSTs) and geographical information systems (GISs) have played a key role for a) optimizing the number of sampling sites for the development of representative sampling networks based on the detection of spatial soil differences and crop growth anomalies [26][27][28], b) the identification of the possible origin of such differences/anomalies through in situ observations/samplings, and c) the design of faster and more targeted interventions.On the other hand, it is very important to stress the high cost of such actions in large agricultural areas and the need for additional indirect methods in order to further reduce their cost through even more targeted soil samplings and soil management interventions.
The aim of this study is to develop an integrated approach to soil quality and fertility assessment in high-yielding rice agro-ecosystems threatened due to overexploitation of soil resources by intensive agriculture.Considering the urgent need for soil monitoring, quality evaluation, and management in rice agro-ecosystems, this approach could contribute to decision-making by farmers and management authorities.
A variety of methodologies, including remote sensing, geostatistics, and multiregression analysis were applied as nodes of an integrated approach.The overall approach was applied to extended pilot fields assuming that soils of similar general properties present a similar nutritional status in large agro-ecosystems where common cropping management practices are followed for a long-term period.The robustness of the aforementioned assumption is indirectly evaluated through measures of statistical significance that derive from the proposed methodologies.The proposed approach is presented in this study using as a case study the large rice agro-ecosystem of Thessaloniki plain in Greece.

Study Area
The study area is the agricultural plain downstream of Axios river in Greece (Figure 2).The mean annual precipitation and temperature in the study area are 440 mm and 15.8 • C, respectively, according to the Worldclim database [29,30], while the mean annual reference evapotranspiration according to the ASCE-standardized method for short grass [31] is 1040 mm [32].The climate of the study area is classified as Cfa/BSk that borders to Csa according to the revised Köppen Geiger classification system [33].The specific area is irrigated by the dam of Axios river.The water is supplied to the fields through an open and extensive collective irrigation network, while drainage is obtained through a drainage network that discharges the water to the sea of Thermaikos Gulf using pump stations.according to the Worldclim database [29,30], while the mean annual reference evapotranspiration according to the ASCE-standardized method for short grass [31] is 1040 mm [32].The climate of the study area is classified as Cfa/BSk that borders to Csa according to the revised Köppen Geiger classification system [33].The specific area is irrigated by the dam of Axios river.The water is supplied to the fields through an open and extensive collective irrigation network, while drainage is obtained through a drainage network that discharges the water to the sea of Thermaikos Gulf using pump stations.The main cultivated crop is rice, where, in most of the cases, a three-year rotation program (maize, cotton, fodder) is applied, while monoculture is usually performed in cases of high-salinity soils.The rice fields are conventionally cultivated with mold board ploughing at November, land leveling at the end of April, and harrowing for seedbed.The irrigation program in the area for rice consists of continuous submergence with short-term interruptions [34,35].
The general nitrogen fertilization practice in these rice fields consists of a) a basal dose of 160 kg N ha −1 in the form of urea, which is provided by broadcasting 500 kg ha −1 of a 32-5-5 chemical fertilizer, and b) a dose of 110 kg N ha −1 as urea at topdressing.This practice has been adopted by the farmers because the new dwarf cultivars used in this area are more nitrogen-demanding and highly productive, reaching yields of 9.5 Mg ha −1 on average (ranging between 8 and 12 Mg ha −1 ) and covering an area of ~19,645 ha based on records of the Hellenic Statistical Authority [36].The specific area constitutes a representative example of a rice agro-ecosystem where rice cultivation has been applied for more than 60 years (the first experimental rice fields appeared in the area in 1949 [37]).The main cultivated crop is rice, where, in most of the cases, a three-year rotation program (maize, cotton, fodder) is applied, while monoculture is usually performed in cases of high-salinity soils.The rice fields are conventionally cultivated with mold board ploughing at November, land leveling at the end of April, and harrowing for seedbed.The irrigation program in the area for rice consists of continuous submergence with short-term interruptions [34,35].
The general nitrogen fertilization practice in these rice fields consists of a) a basal dose of 160 kgN ha −1 in the form of urea, which is provided by broadcasting 500 kg ha −1 of a 32-5-5 chemical fertilizer, and b) a dose of 110 kg N ha −1 as urea at topdressing.This practice has been adopted by the farmers because the new dwarf cultivars used in this area are more nitrogen-demanding and highly productive, reaching yields of 9.5 Mg ha −1 on average (ranging between 8 and 12 Mg ha −1 ) and covering an area of ~19,645 ha based on records of the Hellenic Statistical Authority [36].The specific area constitutes a representative example of a rice agro-ecosystem where rice cultivation has been applied for more than 60 years (the first experimental rice fields appeared in the area in 1949 [37]).
From the overall study area, dispersed pilot rice fields with a total coverage of 753 ha (about ~4% of the area of rice fields) were used in order to perform soil samplings.The pilot rice fields belong to 15 different farmers (Figure 2).

Image Analysis
Image analysis was used in order to detect crop variability over critical stages of rice crop growth in the pilot rice fields, based on the rationale that imagery would reveal systematic soil variation under steadily uniform farming practices followed by the farmers in the past, thus dividing the pilot fields into zones.By default, zonation of fields constitutes a precision agriculture approach; however, consultation, application, and assessments (which are currently underway in the area) were outside the scope of this research.
Sentinel-2 is ready-to-use, medium-high-resolution imagery that is freely available from the European Space Agency (ESA).Its MultiSpectral Instrument (MSI) carries 13 spectral bands ranging from visible and near-infrared to short-wave infrared.The spatial resolution varies from 10 to 60 m, depending on the spectral band, with a 290 km field of view.This ready-to-use product with a unique combination of a high spatial resolution, a wide field of view, and broad spectral coverage in an atmospherically corrected reflectance mode, has opened a new window in operational land use monitoring and mapping.Sentinel-2 satellite imagery of the study area was acquired at three rice growth stages during the 2017 growing season (one on 28 June 2017, one on 28 July 2017, and one on 27 August 2017).The three corrective bands (B1, B9, and B10) were removed as they are not useful for vegetation studies, while the 20-m bands (B5, B6, B7, B8, B8a, B11, and B12) were resampled to 10-m to fit with the rest of the 10-m resolution bands (B2, B3, and B4).Finally, a total of 30 image bands at a 10-m pixel size underwent image segmentation, a process of image division into spatially continuous, disjoint, and homogeneous spatial units called image objects [38].
Image segmentation was performed with the FNEA algorithm (Fractal Net Evolution Assessment) embedded in the eCognition© program (by Definiens Developer).With FNEA, image objects are created by grouping single pixels through a stepwise region-growing procedure of minimizing the internal heterogeneity of each object until a user-defined threshold (called the 'scale factor') is reached [39].The scale factor influences the desired object size and shape [40], according to one spectral and two shape criteria (namely, compactness and smoothness).For segmenting the Sentinel-2 multi-temporal composite, the scale factor was set to 40, while the shape criteria were kept at a low weighting level (0.1), as soil parameters are not expected to impact zone shaping.
The resulted objects divided the land into zones (Figure 3), which are considered to correspond mainly to soil variation, whereas a small partition could be attributed to farm management or other reasons.However, very small and long-shaped zones, specifically, polygons with their attributes AREA < 3000 m 2 or SHAPE INDEX > 2, were filtered out; and, occasionally, filtering was further refined visually in accordance with the affordable number of total soil samples per farmer.Note that, where possible, recent multispectral images of naked soil captured by drones were also used (in addition to a satellite time series).Finally, a total number of 537 sampling positions was set based on the respective number of zones, with the minimum, mean ± standard deviation, and maximum area coverage equal to 0.3, 1.4 ± 0.8, and 6.8 ha, respectively.

Soil and Topographic Parameters
A total number of 537 disturbed soil samples was obtained during the period of February-March 2018 from the 0-30 cm depth of the sampling positions indicated from the previous analysis.The sampling positions and their coordinates were identified in the fields using the mobile app SW-MAPS (freeware).
The soil samples were mechanically crumbled, air-dried, sieved with a 2-mm sieve, and prepared for assessing soil texture through particle size analysis with the hydrometer (or Bouyoucos) method [41].The distribution of samples based on the 12 USDA soil texture classes is given in Figure 4.Additional measurements included: percent of organic matter with the wet digestion method [42], percent of CaCO3 with the back titration method [43], pH measured in saturated paste, and electrical conductivity (EC) measured in a soil saturation extract.Nitrate-nitrogen was measured using 2M KCl [44], available phosphorus according to the Olsen method [45], exchangeable Ca, K, and Mg with ammonium acetate at pH = 7.0 [46], Cu, Mn, Fe, and Zn measured using DTPA [47], and extractable boron with the azomethine-H method [48].
Twenty-one environmental parameters were gathered in total (Table 1) taking also into account topographic parameters obtained by global positioning system (GPS), such as latitude, longitude, and altitude.An additional parameter that describes the distance from the shoreline was also estimated for each sampling position in the GIS environment.The parameters were also divided into two groups (Table 1), where Group 1 includes general soil quality and topographic parameters and Group 2 includes nutrient-related parameters.

Soil and Topographic Parameters
A total number of 537 disturbed soil samples was obtained during the period of February-March 2018 from the 0-30 cm depth of the sampling positions indicated from the previous analysis.The sampling positions and their coordinates were identified in the fields using the mobile app SW-MAPS (freeware).
The soil samples were mechanically crumbled, air-dried, sieved with a 2-mm sieve, and prepared for assessing soil texture through particle size analysis with the hydrometer (or Bouyoucos) method [41].The distribution of samples based on the 12 USDA soil texture classes is given in Figure 4.Additional measurements included: percent of organic matter with the wet digestion method [42], percent of CaCO 3 with the back titration method [43], pH measured in saturated paste, and electrical conductivity (EC) measured in a soil saturation extract.Nitrate-nitrogen was measured using 2M KCl [44], available phosphorus according to the Olsen method [45], exchangeable Ca, K, and Mg with ammonium acetate at pH = 7.0 [46], Cu, Mn, Fe, and Zn measured using DTPA [47], and extractable boron with the azomethine-H method [48].
Twenty-one environmental parameters were gathered in total (Table 1) taking also into account topographic parameters obtained by global positioning system (GPS), such as latitude, longitude, and altitude.An additional parameter that describes the distance from the shoreline was also estimated for each sampling position in the GIS environment.The parameters were also divided into two groups (Table 1), where Group 1 includes general soil quality and topographic parameters and Group 2 includes nutrient-related parameters.  1 Soil organic matter (%) = 2 × Soil organic carbon (%) [49]; 2 The exchangeable Ca was extremely high (above 2000 ppm) in all sampling sites and for this reason was excluded from all types of analysis.  Soil organic matter (%) = 2 × Soil organic carbon (%) [49]; 2 The exchangeable Ca was extremely high (above 2000 ppm) in all sampling sites and for this reason was excluded from all types of analysis.

Soil Quality and Nutrient Deficiency Thresholds
This section was developed in order to provide critical thresholds of the measured soil parameters, which could indicate a high risk of soil quality problems and nutrient deficiencies.The literature that was used for setting these thresholds (except for the case of soil texture) was mostly that of Senadhira [50] and Dobermann and Fairhurst [51], which specialize on rice based on the international literature, and the work of Koukoulakis [52], which provides general thresholds for crops in Greece.The final adopted thresholds were a combination of the aforementioned works taking into account that the rice fields of the study area are described by moderate/highly calcareous soils with very high yields.For this reason, the higher values of the proposed critical limits for nutrient deficiencies in the aforementioned literature were used in order to be more representative of the specific local conditions.

Soil Texture and Organic Matter
Soil texture is probably the most important element that defines the general inherent capacity of soil to regulate the hydraulic properties (water retention, soil aeration, etc.) and the retention and the availability of nutrients.
According to Aschonitis et al. [53], the mean value of long-term yield of winter wheat in different soils of the Thessaloniki region is strongly associated with clay percentage.The yields showed a bell-shaped response versus the increase of clay, where values below <20% and values >50% were responsible for an ~25% yield decrease in comparison to the maximum observed mean yields (see Figure 2 in [53]).The yield decrease for low clay values is probably connected to a lower cation exchange capacity (a low degree of nutrient retention) and a lower water holding capacity and availability, while the respective decrease for the high clay values is probably associated with low soil aeration, problems in water movement, and low water availability due to the higher degree of moisture retention at lower water potentials, which make it difficult for plants to uptake water.In the case of rice, the high values of clay do not create such problems due to the fact that this crop grows under conditions of soil flooding and mud development.On the other hand, low values of clay can create similar nutritional problems in rice as in wheat due to limited retention and availability of nutrients.
An additional threshold related to soil texture can be developed considering sand.Taking into account the hydraulic properties of the 12 USDA soil texture classes according to the two extensive datasets of Carsel and Parrish [54] and Schaap et al. [55], the saturated hydraulic conductivity K s shows a steep increase in the Sandy Loam (SL), Loamy Sand (LS), and Sand (S) classes (more than double K s compared to the other classes; see the complete tables of the two datasets in Aschonitis and Antonopoulos [56]).These three classes correspond to sand% values above 50%.The steep increase of K s above the specific sand% threshold is expected to cause great nutrient leaching and consequently a high risk of fertilization failure.
The aforementioned thresholds of clay <20% and sand >50% can be used either separately or in combination in order to define a high risk of fertilization losses and yield failure in rice fields.The specific thresholds are also verified by Lal [57], who mentioned that the low fertility of LS and SL soils with clay <20% can be easily improved by increasing their soil organic carbon (SOC) pool.
The critical limit of total SOC content, below which crop yield declines by about 20%, is 1.1% for most soils of the tropics [58] and 2% for soils of the temperate regions [57,59].These values correspond to the range of 2.2-4.0% of organic matter (OM% ≈ 2 × SOC%, [49]).For Greek conditions, enrichment of SOC pools is necessary when OM falls below the threshold of 2% [52].In this study, the threshold of OM < 2% was used as a soil quality threshold where a potentially serious decline in soil quality and yield occurs.An upper threshold for OM could also be considered, since high OM values are associated with N, P, Zn, and Cu deficiencies, or Fe and H 2 S toxicities [50].In this study, an upper critical threshold for OM was not used because the observed OM values in the study area are on average 2-3% with rare large values that do not exceed 8% (Table 1).OM values >4% may temporarily occur when straw or manure is incorporated into the soil.

Soil EC, pH, and CaCO 3
The negative effect of salinity stress on rice yields starts to appear for EC values of the saturation extract above 2 mS cm −1 [50,51], while significant negative effects due to soil pH occur for values >7.5 [50].Senadhira [50] gives also a low threshold of pH < 4-5, but this was not used since values below 7.12 were not observed in this study (Table 1).
The EC and pH of irrigation water also plays significant role, and, when its EC exceeds the threshold of 0.5 mS cm −1 and its pH value is outside the range of 6.5-8.0,crop problems start to appear [51].The threshold in the EC value of irrigation water has been set at a significantly lower value compared to the EC value of the soil saturation extract due to the fact that EC of the flooding water in the rice field will significantly increase due to the high evaporation rates.Problems in the study area due to the quality of irrigation water are not expected, since two monitoring programs during 2006-2007 [35] and during 2011 [60] showed that the EC of irrigation water ranges between 0.4 and 0.6 mS cm −1 and the pH ranges between 6.7 and 8.3, indicating a generally good quality of irrigation water.Moreover, the sodium absorption ratio (SAR) values of the specific water are usually below the critical threshold of 15 during the flooding season [35].
Considering the above, the aforementioned thresholds of EC for the soil saturation extract (2 mS cm −1 ) and soil pH (>7.5) are used in this study to evaluate the soil quality status of the specific rice fields without further considerations related to irrigation water quality.
As regards CaCO 3 , values above 5-10% are considered high [52].Even though liming is an appropriate measure for correcting pH in soils, the high CaCO 3 values in non-acidic and calcareous soils may create Fe, Cu, Zn Mn, or B deficiency problems and decrease of P availability due its immobilization in Ca phosphates [51,61].In this study, the value of 7.5% (the average value of 5-10%) was set as the final critical threshold for denoting CaCO 3 as a likely possible limiting factor of rice yield.

Nitrogen and Phosphorus
In irrigated rice systems, the commonly used soil tests are not reliable for predicting soil N supply under flooded field conditions, and therefore they are usually avoided to provide precise critical levels or ranges of optimum nitrogen concentrations in the soil for a rice crop [51].On the other hand, for the purposes related to precision agriculture and fertilization recommendations, thresholds should be defined in order at least to avoid overuse of nitrogen fertilizers.In the case of cereal crops (mainly maize), a range of 20-30 ppm of N in soil is considered to be sufficient [62].For Greek conditions, the specific range for general purposes is expanded to 20-40 ppm [52] due to the semi-arid climate that is usually linked to small organic nitrogen pools [63,64], which could be used as alternative nitrogen sources through mineralization.In this study, the upper threshold of possible nitrogen deficiency for the case of rice was used (i.e., 40 ppm) for the following reasons, taking also into account the suggestions of Dobermann and Fairhurst [51] about possible causes of N deficiency:

•
The specific fields are continuously flooded with a large amount of nitrogen losses through water losses by surface and subsurface drainage [34] in comparison to other non-flooded cereal crops (e.g., maize) for which the aforementioned concentration ranges are mostly applied.High drainage rates are also applied as a strategy to avoid salinity build up though evapo-concentration [35,60].

•
In the specific fields, there are additional losses, such as gaseous losses of nitrogen (volatilization and denitrification) and nitrogen uptake by algae [65][66][67].

•
The total N uptake of rice for grain yields above 10 Mg ha −1 exceeds 200 kg ha −1 [51], and, considering the aforementioned indirect losses in the specific rice agro-ecosystems and their high yields, a higher critical N threshold of serious nitrogen deficiency is required.
As regards phosphorous (Olsen-P), for lowland rice soils with little or no free CaCO 3 , the crop will show a positive response to P fertilization at high yield levels (>8 t ha −1 ) even when the soil P concentration exceeds the value of 10 ppm.In the case of calcareous soils, the respective critical level is raised and possible phosphorus deficiency can occur for Olsen-P concentrations <25 ppm [51].Due to the high yields and the relatively high CaCO 3 concentrations in the rice fields of the study area, the critical threshold for serious P deficiency was set in this study at 25 ppm of Olsen-P.

K, Ca, and Mg
Rice fields with high yield levels (>8 t ha −1 ) can show a positive response to K fertilization even when the K concentration exceeds the value of 175 ppm [51].Ca deficiency is likely to occur for values <200 ppm.The measured values of exchangeable Ca values in all sampling positions of this study were >2000 ppm (Table 1) and for this reason, this element was neglected from the analysis.Mg deficiency is likely to occur for values <120 ppm [51].Thus, in this study, the following critical thresholds were adopted for K, Ca, and Mg deficiency: <175 ppm for K, <200 ppm for Ca, and <120 ppm for Mg.2.4.5.Fe, Zn, Mn, Cu, and B The critical threshold for the occurrence of Fe deficiency in rice fields is <4-5 ppm [51].For Greek conditions and for general use, the critical threshold for possible Fe deficiency is <4 ppm [52].
The critical threshold for the occurrence of Zn deficiency is <0.8-2.0 ppm.For calcareous soils (pH > 7) with moderate to high organic matter (SOC > 1.5%), the upper threshold (2 ppm) should be used because the deficiency is enhanced by the presence of large amounts of HCO 3 − in the soil solution [51].For Greek conditions and for general use, the critical threshold for possible Zn deficiency is <2 ppm [52].
The critical threshold for the occurrence of Mn deficiency in rice is <12-20 ppm [51].For Greek conditions and for general use, the critical threshold for possible Mn deficiency is <15 ppm [52].
The critical threshold for the occurrence of Cu deficiency in rice is <0.3 ppm [51].For Greek conditions and for general use, the critical threshold for possible Cu deficiency is <0.9 ppm [52].
The critical threshold for the occurrence of B deficiency in rice is <0.1-0.7 ppm [51].For Greek conditions and for general use, the critical threshold for possible B deficiency is <0.5 ppm [52].
Taking into account the high rice yields and/or the extensive presence of moderate calcareous soils in the study area, the larger values of the aforementioned critical thresholds for Fe, Zn, Mn, Cu, and B deficiency were adopted (i.e., for Fe < 5 ppm, for Zn < 2 ppm, for Mn < 20 ppm, for Cu < 0.9 ppm, and for B < 0.7 ppm).

Spatial Patterns of Soil Properties
In order to assess the spatial patterns of soil properties in the study area, the statistic of spatial autocorrelation was applied for each soil parameter measured in the laboratory.The Inversed Euclidean Distance method was applied for the tests, revealing that all measured soil properties were significantly clustered (i.e., similar values were found close to each other).These results indicated that the sampling scheme was dense enough to avoid biases due to possible rapid changes in soil properties' values or unexpected discontinuities.Spatial stationarity [68,69] was also checked for all soil properties using global trend tools and was found to be slightly positive over the study area, as expected.Spatial stationarity is defined as the degree of similarity in the way measured variables are autocorrelated throughout the study area; the most non-stationary properties were found to be soil texture, EC, Zn, and Mn (further details would be outside the scope of this research).Furthermore, the clusters were spatially visualized with Hot Spot Analysis, thus revealing the distribution of high and low values over the study area for every parameter.The Hot Spot Analysis was performed according to the Getis-Ord Gi* statistic.For both analyses, the Spatial Statistics toolsets embedded in ArcGIS Pro (ESRI) software were employed [70,71].

Gradient Analysis
All environmental parameters except pH (which is already a log function) were transformed to reduce normality departures (Table 1).Environmental parameters, which are ratios/percentages, were transformed using arcsin(x/100) 0.5 , while the rest were transformed using log(x + 1) [72].The Group 1 variables of Table 1 constitute the "descriptor variables", representing distinct gradients of topography and general soil features; and Group 2 variables constitute the "target variables", representing parameters related to the soil nutritional status.Group 1 variables were used as independent variables and Group 2 variables as dependent variables in order to investigate the degree and the type of nutritional status dependence by general topographic and soil features.Group 1 variables were subjected to a collinearity analysis through the variance inflation factor (VIF).Collinear variables with VIF > 3 were excluded [73] from the following three types of gradient analysis:

•
Multiple gradients analysis: for assessing multiple inter-relations between multiple descriptor variables and multiple target variables.
• Single gradient analysis: for assessing the effect of a single descriptor variable on multiple target variables.

•
Multiple regression: for assessing the effect of multiple descriptor variables on a single target variable.
2.6.1.Multiple Gradients Analysis: Ordination Methods and Variance Partitioning Detrended Correspondence Analysis (DCA) based only on the target variables was used to select the appropriate response model for multiple gradient analysis [74,75].The usual models for multiple gradient analyses are Redundancy Analysis (RDA) (a linear method) and Canonical Correspondence Analysis (CCA) (a unimodal method).When the larger gradient in the results from DCA is below 3 then RDA is more appropriate; when it is larger than 4, CCA is used [75].In this study, the preliminary analysis of target variables with DCA revealed that the dominant gradient length was 0.553 (less than 3), and for this reason RDA was selected for multiple gradient analysis.Descriptors' effect was assessed through variance partitioning based on the marginal effect (λ−1) and the conditional effect (λ-A) of each descriptor variable.The marginal effect of a descriptor variable is equal to the eigenvalue of a partial-RDA if the corresponding variable is the only descriptor variable (additionally to the variance explained by co-variables).The conditional effect of a descriptor variable is equal to the additional amount of variance in the assemblages of target variables explained by the corresponding descriptor variable at the time that was included into the model during a selection procedure [74,76].RDA was performed using CANOCO 4.5 (Biometris -Plant Research International, Wageningen, Netherlands), based on target variable correlations and their standardized scores [74].Significant descriptors were identified using CANOCO's forward selection procedure and Monte Carlo permutation test (499 permutations) (a default option in the CANOCO software).

Single Gradient Analysis: TITAN Method
TITAN (Threshold Indicator Taxa Analysis) [77] is used in ecological studies for detecting changes in taxa abundance distributions along a unique environmental gradient and for assessing synchrony among taxa abundance change points as evidence of community thresholds.In this study, target nutrients' concentration is used instead of taxa abundance in TITAN.TITAN can be used a) to distinguish negative and positive responses and tracks of nutrients' concentrations along a unique gradient (i.e., a general soil property or a topographic property), and b) to detect changes in the distribution of nutrients' concentration and to assess synchrony among their change points along the respective gradient.TITAN uses bootstrapping for estimating purity and reliability (parameters of statistical significance) as well as uncertainty of change points related to individual nutrients' concentrations along the gradient.The cut-off value of 95% was used in both purity and reliability criteria for identifying statistical robust responses of target variables versus each gradient.The purity cut-off value defines what is considered a pure response direction.The default (0.95) requires that 95% of the results from bootstrap replicates agree with the observed response direction.The reliability cut-off value defines what is considered a reliable response magnitude.The default (0.95) requires that 95% of the results from bootstrap replicates have an IndVal p-value less than or equal to 0.05, indicating a response magnitude at a given change point location that differs significantly from what we would expect from random permutation [77] (for explanations about the IndVal p-value, see [77,78]).
TITAN was performed with TITAN2 version 2.1 [79] in the R language using 250 random permutations of targets' variables data and creating 500 new bootstrap datasets by resampling the observed data with a replacement (default options in TITAN2).The analysis was made only for the important descriptor variables (gradients) indicated by RDA analysis, and the results are provided as response plots of the pure and reliable (at the 95% level) target variables versus each gradient.The meaning of this type of analysis is to detect if there are specific points (i.e., a value) along the gradient (i.e., a range of values) of a descriptor variable where the nutrient concentrations (target variables) present significant changes (positive or negative).

Multiple Linear Regression
Multiple Linear Regression (MLR) was performed in two steps, which are described in the following two paragraphs.
In the first step, a typical MLR for each target variable was performed using forward selection of descriptor variables.In this way, the statistically significant descriptor variables were selected in the final models.The final MLR models for each target variable were used for repeating MLR with bootstrapping (Boot-MLR).Bootstrapping is used to generate a large number of new datasets by randomly sampling data with a replacement [80] and it is considered to be among the most robust methods for assessing the variability and robustness of regression coefficients.It is also an alternative procedure that exploits the whole dataset, avoiding data split into calibration and validation groups for verifying models' performance.Boot-MLRs were performed by applying the "nls.lm"function of the {minpack.lm}package [81] in R software.The "nls.lm" function uses the Levenberg-Marquardt least-squares algorithm.The Boot-MLR procedure was applied for 10,000 iterations that led to a respective number of solutions for the regression coefficients.The range of bootstrap estimations of the regression coefficients was defined by the 95% confidence interval, which was estimated based on the probability distribution of their 10,000 estimations.This method was applied in order to estimate the values of the highest posterior density distribution (HPDD) that indicates the 2.5% and 97.5% thresholds (HPDD thresholds), which contain the central 95% of the HPDD.The 95% HPDD probability interval was computed using the "p.interval" of the {LaplacesDemon} package [82] in R software.This analysis was performed in order to investigate the robustness of the coefficients, and it was assumed that if the coefficients preserve their positive or negative signal within the 95% HPDD interval, they are considered robust.
In the second step, the selected descriptors of the typical MLR models for each target variable were used for repeating Ridge MLR with standardized values.In this way, the values of estimated regression coefficients could provide a comparison of the relative effect of each descriptor variable on a target variable.Ridge MLR analyses were performed using StatGraphics Centurion XV software (StatPoint Inc., The Plains, VA, USA).

Spatial Patterns of Soil Parameters and Quality/Nutritional Status of Rice Soils
The percentage of the corresponding area of pilot fields that violate the quality thresholds of Section 2.4 were calculated individually for each parameter and they are given in Table 2.According to Table 2, individual deficiencies of exchangeable Ca, Mg, and micronutrients, such as Fe and Cu, do not exist in almost ~100% of the pilot total area.
The most important problems in terms of percent coverage (greater than 50% of the pilot fields' area) are the large individual occurrence of higher values of pH and CaCO 3 above the critical thresholds and the large individual occurrence of nitrate-nitrogen, phosphorus, zinc, manganese, and boron below the critical thresholds.The most important combined violations of critical thresholds in terms of percent area coverage are the following:

•
The regions of simultaneous deficiency in Zn, Mn, and B cover 62.5% of the total area of pilot fields.

•
The regions of simultaneous pH > 7.5 and CaCO 3 > 7.5% cover 39.7% of the total area of pilot fields.

•
The regions of simultaneous deficiency in N, P, and K cover 22.5% of the total area of pilot fields.
The aforementioned observations and other spatial patterns of each parameter of Table 2 were visualized through hot/cold spot analysis using the Getis-Ord Gi* statistic, which provides significant clusters of high values and low values for each parameter (Figure 5).

Multiple Interrelations between Descriptor Variables and Target Variables: RDA Analysis
The descriptor variables lat, long, alt, S, and Si were excluded through forward selection of RDA analysis in order to form a group of descriptor variables with VIF values and high statistical significance.The finally selected descriptor variables (OM, C, EC, pH, dist, and CO3) all showed VIF

Multiple Interrelations between Descriptor Variables and Target Variables: RDA Analysis
The descriptor variables lat, long, alt, S, and Si were excluded through forward selection of RDA analysis in order to form a group of descriptor variables with VIF values and high statistical significance.The finally selected descriptor variables (OM, C, EC, pH, dist, and CO3) all showed VIF < 2.7 and statistical significance at the p < 0.01 level.The importance of descriptor variables based on their marginal (λ-1) and conditional (λ-A) effects is visualized in Figure 6a, while the triplot between the descriptor variables or "gradients", the target variables, and the sampling sites provided by RDA is provided in Figure 6b.Overall, the proportion of variance explained by the selected descriptor variables (OM, C, EC, pH, dist, and CO3), expressed as the sum of all canonical eigenvalues, was 49.4%.The majority of variance in the targets-descriptors relation is mostly explained by Axis 1 with 71.4%, while Axis 2's contribution is restricted to 25.2% (the remaining 3.4% is explained by Axes 3 and 4, which are not visualized in Figure 6a).
Taking into account Figure 6a, the importance of descriptor variables on NO3N, P, K, Fe, Mg, Mn, Cu, Zn, and B variation based on their marginal and conditional effect follows the order OM > C > EC > pH > dist > CO3.Taking into account Figure 6b, the most important observations are the following:

•
C and OM values are negatively correlated with dist and NO3N values and positively correlated with P, K, Fe, Mg, Mn, Cu, Zn, and B. Thus, soils closer to the coastline are more likely to have higher C and higher OM, while they also have a lower probability to present P, K, Fe, Mg, Mn, Cu, Zn, and B deficiency.On the other hand, the soils with lower OM and C values that are far from the coastline are more likely to show higher values of NO3N and a higher probability to present P, K, Fe, Mg, Mn, Cu, Zn, and B deficiency.
EC is negatively correlated with pH and positively correlated with P. Thus, soils of higher EC are more likely to show higher concentrations of P and lower pH values.
The aforementioned elements can also be spatially located by combining the maps of Figure 5.
Agriculture 2019, 9, x FOR PEER REVIEW 15 of 26 from the coastline are more likely to show higher values of NO3N and a higher probability to present P, K, Fe, Mg, Mn, Cu, Zn, and B deficiency.• EC is negatively correlated with pH and positively correlated with P. Thus, soils of higher EC are more likely to show higher concentrations of P and lower pH values.The aforementioned elements can also be spatially located by combining the maps of Figure 5.

Single Gradient Analysis: Association of a Single Descriptor Variable with Multiple Target Variables
The most important descriptor variables from the RDA (OM, C, EC, pH, dist, and CO3) (Figure 6a) were further used in order to detect their associations with the target variables using the TITAN method.The association of each one single descriptor variable versus multiple pure target variables according to TITAN is provided in Figure 7.
For each descriptor variable in Figure 7, only pure indicator target variables are plotted in increasing order with respect to their observed environmental change point.Those showing a negative response to the gradient are on the left side of the plot (black-filled circles with black lines), while those with a positive response are on the right side (white-filled circles with dashed lines).Circles are sized in proportion to Z scores defining the most significant change points, while the horizontal lines overlapping each circle represent the 5th and 95th percentiles of change points

Single Gradient Analysis: Association of a Single Descriptor Variable with Multiple Target Variables
The most important descriptor variables from the RDA (OM, C, EC, pH, dist, and CO3) (Figure 6a) were further used in order to detect their associations with the target variables using the TITAN method.The association of each one single descriptor variable versus multiple pure target variables according to TITAN is provided in Figure 7.
For each descriptor variable in Figure 7, only pure indicator target variables are plotted in increasing order with respect to their observed environmental change point.Those showing a negative response to the gradient are on the left side of the plot (black-filled circles with black lines), while those with a positive response are on the right side (white-filled circles with dashed lines).Circles are sized in proportion to Z scores defining the most significant change points, while the horizontal lines overlapping each circle represent the 5th and 95th percentiles of change points among 500 bootstrap replicates.Figure 7a shows that P, K, Mg, Mn, Zn, Fe, Cu, and B present a positive response to the increase of organic matter with significant positive change points (white-filled circles on the dashed lines) that appear for OM >1.9% (the minimum value that corresponds to Mg) and up to 2.8% (the maximum value that corresponds to Mn).Similarly, Figure 7b shows that all of the aforementioned nutrients show a positive response to clay increase with significant positive change points that appear for C >22% (the minimum value that corresponds to Cu) and up to 40% (the maximum value that corresponds to B).In both Figures 7a,b, NO3N shows a negative response to organic matter and clay increase with significant negative change points that occur for OM = 2.6% and C = 22%.
In Figure 7c, observed is a positive response of P, K, Mg, Mn, Zn, Fe, Cu, and B with the increase of electrical conductivity with significant positive change points that appear for EC >0.7 mS cm −1 (the minimum value that corresponds to Mg) and up to 2 mS cm −1 (the maximum value that corresponds to P).
Figure 7d shows that P, K, Mg, Mn, Zn, Fe, Cu, and B present a negative response to the increase of pH with significant negative change points (black-filled circles on the black lines) that appear for Figure 7a shows that P, K, Mg, Mn, Zn, Fe, Cu, and B present a positive response to the increase of organic matter with significant positive change points (white-filled circles on the dashed lines) that appear for OM > 1.9% (the minimum value that corresponds to Mg) and up to 2.8% (the maximum value that corresponds to Mn).Similarly, Figure 7b shows that all of the aforementioned nutrients show a positive response to clay increase with significant positive change points that appear for C > 22% (the minimum value that corresponds to Cu) and up to 40% (the maximum value that corresponds to B).In both Figure 7a,b, NO3N shows a negative response to organic matter and clay increase with significant negative change points that occur for OM = 2.6% and C = 22%.
In Figure 7c, observed is a positive response of P, K, Mg, Mn, Zn, Fe, Cu, and B with the increase of electrical conductivity with significant positive change points that appear for EC > 0.7 mS cm −1 (the value that corresponds to Mg) and up to 2 mS cm −1 (the maximum value that corresponds to P).
Figure 7d shows that P, K, Mg, Mn, Zn, Fe, Cu, and B present a negative response to the increase of pH with significant negative change points (black-filled circles on the black lines) that appear for pH > 7.8 (the minimum value that corresponds to Zn) and up to 8.0 (the maximum value that corresponds to Cu).
The association of the distance from the shoreline dist with K, Mg, Mn, Zn, Fe, Cu, and B introduces the aspect of setting topographic thresholds regarding nutrient deficiencies.Figure 7e shows that all of the aforementioned nutrients present a negative response to the increase of dist with significant negative change points that appear in the dist range between 1 and 2.5 km from the shoreline.P did not have a pure association with dist.On the other hand, NO3N shows a positive response to dist increase with a significant change point that occurs for dist = 4.0 km.
Figure 7f shows that K, Mg, Mn, Zn, Cu, and B present a positive response to the increase of CO3 with significant positive change points that appear for CO3 > 6.5% (the minimum value that corresponds to Cu) and up to 10.3% (the maximum value that corresponds to Zn).On the other hand, P and NO3N show a negative response to CO3 increase with significant negative change points that occur for CO3 = 7% and CO3 = 9.2% for P and NO3N, respectively.

Regression Models
In the first step, a typical MLR analysis for each individual target variable (NO3N, P, K, Mg, Fe, Zn, Mn, Cu, and B) using forward selection of descriptor variables (dist, C, OM, CO3, pH, and EC) was performed in order to select the statistically significant descriptor variables for developing the final form of MLR models.The final form of the typical MLR models for each target variable was used for repeating MLR with bootstrapping (Boot-MLR).The results of Boot-MLR for each target variable are given in Table 3.
Based on the results of Table 3, the models were all significant (according to the p-value).The robustness of the Boot-MLR models of each target variable was ranked based on the adjusted R 2 for the respective degrees of freedom (R 2 adj.df ) as follows K > Mg > Cu > P > Mn > Fe > Zn > B > NO3N.In the second step, the selected descriptors of the typical MLR models for each target variable were used for repeating Ridge-MLR with standardized values.The values of the estimated regression Ridge-MLR coefficients are given in Figure 8.
The models, but also the coefficients, of Boot-MLR (Table 3) and Ridge-MLR (Figure 8) are all statistically significant; however, the models' majority presents a medium-to-low predictive accuracy according to R 2 adj.df .This indicates that the models can mostly be used to predict the possible trends and the possible negative or positive effect of a descriptor variable rather than to be directly used for predicting precise values of target variables.Figure 8 provides an alternative perspective of the descriptor variables effect, which is in accordance with the RDA (Figure 6b) and TITAN (Figure 7) results.The additional information that is given by Figure 8, compared to RDA and TITAN, provides a comparison between the relative effects of the descriptor variables on each target variable.
In the case of using the MLR models for predicting the value and not the trend of a target variable, an R 2 adj.df < 0.5 suggests a high risk of predictive accuracy.The MLR models that show R 2 adj.df ≥ 0.5 are those of the K, Mg, Cu, and P parameters (Table 3).Mg and Cu did not show deficiency in ~100% of the pilot total area (Table 2).For this reason, special attention should be paid to the MLR models of K and P, which are macronutrients and constitute basic elements in N-P-K mixtures of fertilizers.

3.
Results of multiple linear regression with bootstrapping (Boot-MLR) models for each individual target variable (NO3N, P, K, Mg, Fe, Zn, Mn, Cu, and B) after forward selection of descriptor variables (dist, C, OM, CO3, pH, and EC) using typical MLR (the MLR models are based on transformed variables according to the transformations of Table

The Role of Each Type of Analysis in Developing Efficient Integrated Interventions for Improving Soil Quality and Productivity
This study provided different types of analysis that can be used to develop effective integrated management schemes.The role of each analysis is described in the following paragraphs.
The critical threshold analysis and the Hot-Cold Spot analysis that were presented in Section 3.1 indicated the most important soil quality and nutrient deficiency problems in space and in terms

The Role of Each of Analysis in Developing Efficient Integrated Interventions for Improving Soil Quality and Productivity
This study provided different types of analysis that can be used to develop effective integrated management schemes.The role of each analysis is described in the following paragraphs.
The critical threshold analysis and the Hot-Cold Spot analysis that were presented in Section 3.1 indicated the most important soil quality and nutrient deficiency problems in space and in terms of area coverage.This element can be used to set priorities for intervention in order to improve the soil quality conditions of the rice fields based on collective actions that could improve drastically the mean productivity of the total study area.For example, 96.6% of the total area showed pH > 7.5.Moreover, 62.5% of the total area showed simultaneous deficiency in Zn, Mn, and B; thus, a combined intervention to improve the aforementioned elements would show an important effect in the mean productivity of the total study area.The areas of intervention can be located by the Hot Spots of pH and the Cold Spots of Zn, Mn, and B given in Figure 5.
The RDA that was presented in Section 3.2 provided a) a ranking of the general soil parameters and topographic elements "descriptor variables" and b) the gradient trends, based on their multiple inter-relations with the nutrient-related target variables.Figure 6a indicates the importance of descriptor variables to describe the spatial variation of target variables.For example, the organic matter is the most important gradient, which means that a collective intervention to improve this parameter to optimum levels would show an important effect in the mean productivity of the total study area due to its positive correlation with the majority of nutrients (except nitrogen, which is anyway provided by the fertilizers).From Figure 5, it is observed that the Cold Spots of OM show important overlapping with the Cold Spots of K, Mg, Zn, Mn, Cu, and B, indicating that the improvement of the organic matter conditions at these Cold Spots could also improve the concentrations of the aforementioned elements and the general soil quality and fertility [83][84][85].
TITAN, which was presented in Section 3.3, provides threshold values of the general soil parameters and topographic elements "descriptor variables", which are significant change points of steep increase or decrease in nutrients' availability.Since organic matter was the most important parameter according to RDA (Figure 6a), an intervention to the Cold Spots of organic matter in order to increase its level to around 2−3% can simultaneously achieve higher rates of availability in P, K, Mg, Mn, Zn, Fe, Cu, and B (Figure 7a).
The regression analysis that was presented in Section 3.4, and more specifically the MLRs, can detect the individual trends of each nutrient-related target variable using the descriptor variables.Moreover, if an MLR shows an adequate predictive accuracy, specific values of the target nutrient-related variable can also be assessed.For example, the MLR models for P and K macronutrients, which showed an adequate performance, can be used to make specific estimations for their concentrations based on data of general soil properties that exist not only for the pilot fields but also for the whole area, since detailed soil surveys and maps [86][87][88], but also soil samples of individual farmers, are analyzed every year in the laboratories of the Soil and Water Institute of the Hellenic Agricultural Organization "Demeter" in Thessaloniki (Greece).Such generalized estimations can be taken into account in order to assess the spatial proportions and the degree of K and P deficiencies in order to propose optimum recipes of N-P-K mixtures including Zn, Mn, and B micronutrients to local fertilizer industries and suppliers.

Interesting Observations Made in the Context of This Case Study
The results of this study revealed significant information about the inter-relations between soil parameters, which raised a lot of questions in the context of soil management.
For example, it was really interesting that nitrate-nitrogen was negatively correlated with clay and organic matter (Figures 6b and 7a,b).Looking carefully the data of sampling positions with the higher nitrate concentrations, it was observed that these soils, apart from low clay and organic matter, had the highest sand percentage and belonged to the Loamy Sand and Sandy Loam soil texture classes.Two possible explanations about this observation are: a) the incorporated fresh organic material (roots, parts of straw) in these light soils from the previous growing season was mineralized faster due to the higher aeration and the low retention of soil water, b) a significant amount of nitrates was leached due to higher hydraulic conductivities below the main root zone during the previous growing season and re-appeared on the surface due to soil inversion/blending during the application of deep mold board ploughing.
Another interesting observation was a positive response of P, K, Mg, Mn, Zn, Fe, Cu, and B concentrations with the increase of electrical conductivity (Figure 7c).Based on personal communication with the farmers, the fields of high salinity usually give 10−20% lower yields in comparison to the fields with no salinity problems.Considering this, it can be inferred that the higher concentrations of the aforementioned elements in the soils of higher EC may be due to the lower degree of uptake by the plants, which are affected by the high osmotic stress.A similar explanation can also be given for the effects of CaCO 3 , where there is a positive response of K, Mg, Mn, Zn, Fe, Cu, and B concentrations (except P) for CaCO 3 > 6% (Figure 7e).In this case, there is probably a lower degree of their uptake by the plants, not due to osmotic stress, but because there is a high disturbance in the ratios among Mg:Ca, Mg:K, Ca:K, and (Ca + Mg)/K that significantly affect crop growth by creating unbalanced nutrient concentrations in the soil solution [51,89], but also because the high CaCO 3 leads to P immobilization in Ca phosphates, restricting phosphorus availability, which is also verified in Figure 7e.

Uncertainties Regarding Quality/Deficiency Thresholds
In Section 2.4, an effort was made to define quality/deficiency thresholds in soil parameters taking into account the literature that is more compatible with the features of the study area (e.g., high yields and moderate calcareous soil conditions).Regarding the yields, it is also important to note that small pilot experiments in lysimeters for estimating rice production potential and water consumption under the local climatic conditions showed that yields can even exceed 15 Mg ha −1 with maximum leaf area index values between 10 and 12 [90,91].These elements should be considered when farmers set high target yields, since the higher yields require higher amounts of nutrients and water.For this reason, there is large uncertainty in the case of nutrient deficiency thresholds, when measurements of nutrients concentrations in the produced biomass are not performed.Such measurements are really important in order to estimate the exact nutrient amount that is extracted per unit of biomass production.Hydro-bio-chemical modeling approaches and in situ measurements are also required in order to estimate the water requirements, and the inflows and outflows of nutrients through other processes (e.g., inflows through irrigation and rainfall, outflows through percolation, runoff, denitrification, ammonia volatilization, algae uptake, etc.) [65][66][67]92].Thus, all the aforementioned information should be combined in order to define deficiency thresholds that take into account the targeted/expected yields and the expected amount of nutrient losses through processes not related to plant uptake.The thresholds that were provided in this study are not for general use and satisfy only a) evaluation purposes for the conditions of the study area and b) the presentation of the proposed combination of methodologies as an integrated approach that can support the development of soil management strategies in large agro-ecosystems.For other case studies, the proposed thresholds should be adapted considering the respective site characteristics (e.g., average yields, soil attributes related to observed problems, and irrigation and drainage strategies).
Other missing information that could significantly support the proposed approach is the final yields per zone.This information requires advanced threshing machines equipped with precision agriculture devices in order to perform yield measurements per zone.This would allow for the inclusion of the term of productivity as an additional target variable in the proposed modeling/statistical approaches of this study and to make more robust adjustments in the quality/deficiency thresholds.
Finally, it is important to address that parameters may have more than one critical threshold.For example, high values of some nutrients can lead to toxicities.In this study, such thresholds where not used due to their generally low concentrations.On the other hand, there are some parameters that require more attention because a critical threshold that defines crop tolerance limits can hide significant information related to other problems.For example, the pH critical threshold of 7.5 revealed that the rice of the study area grows under non-optimal pH conditions (Table 2) but at the same time hindered the existence of more serious problems.For example, 40.5% coverage of the pilot fields showed pH ≥ 8, while one pilot zone of ~1.3 ha showed pH > 8.5.These observations about pH could suggest a trend towards sodic conditions in a significant part of the study area.For a general case, the aforementioned element should also be considered by including additional parameters related to sodicity, such as Na and Exchangeable Sodium Percentage (ESP), in the proposed methodology in order to design relevant management practices.Of course, this problem in this study area is already known and has been studied through various soil surveys [86][87][88] and through soil analyses of individual farmers.For this reason, a detailed analysis was not performed regarding sodicity parameters.The specific area had extremely high salinity and sodicity problems even before the first appearance of rice fields due to the existence of the sea from the ancient times, which was gradually converted to lakes and wetlands due to depositions from rivers, such as Axios, Gallikos, and Aliakmon.Finally, the area was converted to agricultural land through extensive reclamation works that started at the beginning of the 19th century and continued to the mid-1960s, providing the current form of the study area [93].Rice was introduced in 1949 as the best solution not only due to its high productivity in the specific soils but also as a means to reduce salinity through salt leaching by flooding in combination with additional measures, such as gypsum application for reducing sodicity [37].

Conclusions
This study proposed an integrated approach of combined methodologies for assessing the soil quality and nutritional status of large rice agro-ecosystems where common cropping management practices have been followed for a long time period.The application of the proposed approach on the rice agro-ecosystem of Thessaloniki plain was able a) to identify the most important individual and combined soil quality and nutrient deficiency problems in space and in terms of area coverage, b) to rank the importance of general soil properties, such as organic matter, soil texture, pH, salinity, and CaCO 3 , or non-soil parameters (e.g., distance from the coastline) as descriptor gradients for assessing the overall or the individual trends of nutrient-related target variables, and c) to identify threshold values of the descriptor gradients where the nutrient-related target variables show a significant negative or positive change.Based on the results of the aforementioned procedures, the minimum proposed actions that could significantly improve the soil nutritional status of the specific agro-ecosystem and consequently its respective production are (a) the development of targeted interventions to increase organic matter and Zn, Mn, and B micronutrients, and (b) to assess the spatial proportions and the degree of K and P deficiencies in order to propose optimum recipes of N-P-K mixtures including Zn, Mn, and B micronutrients to local fertilizer industries and suppliers.The proposed approach was applied to the specific case study considering some site specificities (e.g., high rice yields, calcareous soils) when setting the critical thresholds in soil parameters.Thus, adaptations in the critical thresholds are required before any application to other study sites.

Funding:
The soil samplings and soil analyses were funded by the farmers.

Figure 2 .
Figure 2. Pilot rice fields of the study area (colored in yellow).

Figure 2 .
Figure 2. Pilot rice fields of the study area (colored in yellow).

Figure 3 .
Figure 3. Delineation of detected zones with image segmentation of satellite image time series and their corresponding soil samples in a group of rice fields in the central sector of the study area.

Figure 3 .
Figure 3. Delineation of detected zones with image segmentation of satellite image time series and their corresponding soil samples in a group of rice fields in the central sector of the study area.

Figure 4 .
Figure 4. Distribution of the 537 samples based on the 12 USDA soil texture classes.

Figure 4 .
Figure 4. Distribution of the 537 samples based on the 12 USDA soil texture classes.

Figure 5 .
Figure 5. Hot Spots (significant clusters of high values) and Cold Spots (significant clusters of low values) according to the Getis-Ord Gi* statistic for each parameter of Table 2. Larger z-scores indicate more intense clustering of high values (hot spot); smaller z-scores indicate more intense clustering of low values (cold spot).

Figure 5 .
Figure 5. Hot Spots (significant clusters of high values) and Cold Spots (significant clusters of low values) according to the Getis-Ord Gi* statistic for each parameter of Table 2. Larger z-scores indicate more intense clustering of high values (hot spot); smaller z-scores indicate more intense clustering of low values (cold spot).

Figure 6 .
Figure 6.(a) Marginal and conditional effects of descriptor variables according to Redundancy Analysis (RDA), (b) a triplot of RDA considering the two major axes.

Figure 6 .
Figure 6.(a) Marginal and conditional effects of descriptor variables according to Redundancy Analysis (RDA), (b) a triplot of RDA considering the two major axes.

Figure 8 .
Figure 8.The coefficients value from Ridge-MLR for each individual nutrient-related target variable.

Table 1 .
Statistics of environmental parameters obtained from sampling positions of pilot fields.

Table 1 .
Statistics of environmental parameters obtained from sampling positions of pilot fields.

Table 2 .
Percentage of the total area of pilot fields that violates the soil quality/nutrient thresholds defined in Section 2.4.
* % total area was calculated based on the area of each zone that corresponds to a specific sampling point.Agriculture 2019, 9, x FOR PEER REVIEW 14 of 26 1, * p < 0.05, ** p < 0.01).The coefficients value from Ridge-MLR for each individual nutrient-related target variable.