**Sustainable Agriculture and Soil Conservation**

Editors

**Mirko Castellini Concetta Eliana Gattullo Anna Maria Stellacci Mariangela Diacono**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin


*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Applied Sciences* (ISSN 2076-3417) (available at: www.mdpi.com/journal/applsci/special issues/ Sustainable Agriculture Soil Conservation).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-2477-1 (Hbk) ISBN 978-3-0365-2476-4 (PDF)**

© 2021 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**



### *Editorial* **Sustainable Agriculture and Soil Conservation**

**Mirko Castellini 1,\* , Mariangela Diacono <sup>1</sup> , Concetta Eliana Gattullo <sup>2</sup> and Anna Maria Stellacci <sup>2</sup>**


**Abstract:** Soil degradation is one of the most topical environmental threats. A number of processes causing soil degradation, specifically erosion, compaction, salinization, pollution, and loss of both organic matter and soil biodiversity, are also strictly connected to agricultural activity and its intensification. The development and adoption of sustainable agronomic practices able to preserve and enhance the physical, chemical, and biological properties of soils and improve agroecosystem functions is a challenge for both scientists and farmers. This Special Issue collects 12 original contributions addressing the state of the art of sustainable agriculture and soil conservation. The papers cover a wide range of topics, including organic agriculture, soil amendment and soil organic carbon (SOC) management, the impact of SOC on soil water repellency, the effects of soil tillage on the quantity of SOC associated with several fractions of soil particles and depth, and SOC prediction, using visible and near-infrared spectra and multivariate modeling. Moreover, the effects of some soil contaminants (e.g., crude oil, tungsten, copper, and polycyclic aromatic hydrocarbons) are discussed or reviewed in light of the recent literature. The collection of the manuscripts presented in this Special Issue provides a relevant knowledge contribution for improving our understanding on sustainable agriculture and soil conservation, thus stimulating new views on this main topic.

**Citation:** Castellini, M.; Diacono, M.; Gattullo, C.E.; Stellacci, A.M. Sustainable Agriculture and Soil Conservation. *Appl. Sci.* **2021**, *11*, 4146. https://doi.org/10.3390/ app11094146

Received: 25 April 2021 Accepted: 27 April 2021 Published: 1 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Keywords:** organic agriculture; soil amendment; soil organic carbon; soil management; soil water repellency; soil contamination; soil remediation; sustainable fruit growing; water conservation practices; multivariate statistical models for SOC prediction

#### **1. Introduction**

The ongoing global climate changes and the ever-increasing world population impose the adoption of agroecological or eco-sustainable production techniques to increase the resilience of the most fragile agrosystems [1]. For this reason, from many directions, a multidisciplinary scientific paradigm is increasingly suggested to better approach the complex problems concerning soil protection and sustainable agriculture [1–3].

The development and adoption of sustainable agronomic practices able to preserve and enhance the physical, chemical and biological properties of soils and improve agroecosystem functions is a challenge for both scientists and farmers [4]. In recent decades, much has been done because, on one hand, scientists have improved and simplified the measurement methods in specific research sectors, and on the other hand, the process of transferring research outputs to farmers has been enhanced. The latter has been considerably extended thanks to research funding granted by the National and Supranational Political Institution, which increasingly involves farmers and stakeholders. From a scientific point of view, much progress has been achieved in the development of sustainable practices based on minimal soil disturbance, the use of cover crops, organic mulching, crop rotations, water and nutrient conservation, recycling crop residues and livestock manure for soil amendment, precision agriculture and so on, and many reviews have summarized the state of the art and traced the future trajectories (among others, refs. [5–9]). However, much remains to be explored, and Special Issues stimulating the scientific community on these main topics are desirable.

The main goal of this Special Issue (SI) is to collect and present recent research on sustainable agriculture and soil conservation. A synthesis of the main results is reported in the following section.

#### **2. Overview of This Special Issue**

This SI collects 12 original contributions focused on soil use and management, soil conservation, as well as on the impact of some pollutants on the soil. Specifically, the SI collects nine research papers [10–18], two reviews [19,20] and one methodological work [21]. Eight of the twelve manuscripts considered the effects of soil use and management (topic I) [10–16,21], while the remaining four [17–20] evaluated the impact of some soil contaminants (e.g., crude oil, tungsten, copper and polycyclic aromatic hydrocarbons) (topic II). This grouping was used to detail the contents and for the reader's usefulness.

Regarding the aforementioned nine investigations (i.e., research papers), seven were carried out in the field [10–16], while the others accounted for pot [17] or microcosm [18] experiments. From a geographical point of view, five "field open-air" investigations were carried out in Italy (three in the south [12,13,16], one in the center [15] and one in the north [10] of Italy), while the remaining investigations were in Brazil [14] and China [11].

Topic I comprises eight papers. Pastorelli et al. [10] investigated the short-term effects of digestate on soil properties through a holistic approach, thus investigating the soil's physical, chemical and microbiological properties and their interactions. The digestate effects were evaluated in the two maize growing seasons within a 2-year maize–triticale rotation to feed the biogas plant. The main results showed that the digestate (1) increased the soil's total organic C, total N and K contents, (2) did not affect the soil bulk density while transiently improving the soil aggregate stability, (3) decreased the soil transmission pores (50–500 µm size) proportion as well as increasing the fissures (>500 µm), (4) only transiently affected the soil's microbial community and did not cause soil contamination from Clostridiaceae-related bacteria, (5) significantly impaired seed germination when applied at low dilution ratios and (6) did not improve the crop yield, as it was similar to usual mineral fertilization. In the context of energy crop farming, the authors concluded that the agronomic recycling and usage of digestate from biogas production assured a sustainable crop yield and soil quality. Therefore, digestate was confirmed as a valid resource for sustainable management of the soil's fertility.

Wang et al. [11] investigated the prediction of soil organic matter (SOM) from visible and near-infrared (VNIR) spectra, taking into account the presence of nonlinear relationships and spatial heterogeneity conditions. Specifically, they combined the proposed partial least squares-based multivariate adaptive regression spline (PLS–MARS) method and a regional multivariable associate rule mining and Rank–Kennard–Stone method (MVARC-R-K-S) to construct a nonlinear prediction model to realize local optimality considering spatial heterogeneity. The method, which was applied and validated for a case study in Hubei (China), was able to filter the spatial autocorrelation present in the investigated area. In addition, it showed improvements compared with some available methods in terms of accuracy and robustness, suggesting the reliability of the proposed prediction model.

Losciale et al. [12] compared four soil management strategies at a peach orchard, namely (1) completely tilled (control); (2) mulched with reusable reflective plastic film; (3) mulching with a leguminous cover crop flattened after the peach fruit set; (4) completely tilled and irrigated with the water volumes of the treatment with reflective mulching. The effects on the soil features, water use efficiency (WUE), tree functionality, fruit growth and quality, yield and water productivity were investigated. The main results showed that when using 50% of the regular irrigation, reusable reflective mulching reduced the water loss and soil carbon mineralization while not affecting (and sometimes increasing) the net carbon assimilation, yield and fruit size. Consequently, water productivity was drastically increased. The investigation allowed for the conclusion that a reflective mulching strategy

is promising from the point of view of WUE optimization, especially in hot and dry areas with clay soils and low organic matter contents.

De Mastro et al. [13] investigated the effects of different soil management methods (conventional tillage, CT, minimum tillage, MT, and no-tillage, NT), both with fertilization and without fertilization, on the quantity of SOC at different soil depths and soil particle size fractions. Diffuse reflectance infrared Fourier transform spectroscopy (DRIFT) was also performed to obtain information on both the SOC quality and the mineralogical composition of these soil fractions. The results showed that CT provided the highest amount of the finest fraction, while fertilization improved the microbial community with the increase of soil microaggregates. The coarse fraction was highest in the upper soil layer, while the finest fraction was highest in the deepest one. The greatest OC content was observed in the topsoil layer and in the finest soil fraction. The DRIFT analysis could provide information about the quality of the main minerals present in the different soil size fractions.

Lozano-Baez et al. [14] measured both the unsaturated and saturated soil hydraulic conductivity, determined with simple and low-cost field infiltration methods—a minidisk tension infiltrometer, MDTI [22], and Beerkan [23]—for three land covers, namely (1) pastures, (2) 9-year-old restored forest and (3) remnant forest. The paper, which was relevant because it provided the first measurements of soil water repellency (SWR) in the Brazilian Atlantic Forest, proved that MDTI and Beerkan infiltrations gave complementary information, highlighting increasing hydraulic conductivities, especially at the remnant forest plots, when moving from near-saturated to saturated conditions. Moreover, measuring the hydraulic conductivity with MDTI allowed the estimation of the macroscopic capillary length [24], which is a hydrodynamic soil parameter useful for investigating the impact of soil management on water saving [25]. The authors concluded that the applied approach (MDTI + Beerkan) allowed the design of better estimates of the saturated soil hydraulic conductivity under challenging field conditions, such as soil water repellency (SWR).

Amoriello et al. [15] evaluated the impact of soil amendment with biochar from brewers' spent grain (BSG) on hop growth. The field investigation was carried out in Rieti (central Italy). Three different German cultivars of hop plant were considered, and biochar was added at a 20% level. Biochar's effects on root development, shoots, bine length and crop yield were evaluated using multivariate image analysis and general linear models. The results showed that biochar significantly improved root growth. Overall, no variability in shoots, number of leaves or bine length was observed between the two treatments for all cultivars or among the genotypes considered. Moreover, soil amendment significantly improved the yield (number of cones). The authors concluded that BSG-derived biochar could be useful for improving hop plant growth and cone production because it can supply key nutrients for plant growth and improve the soil's properties.

Montemurro et al. [16] evaluated the production capacity of organic horticultural systems and the ex-post sustainability using a new multi-attribute decision model named DEXi-met. Specifically, they compared three horticultural systems in Metaponto (south Italy): (1) ECO (organic system with full implementation of agroecological strategies, agroecological services crops (ASC), strip cultivation and organic amendment), (2) GM (organic system with the introduction of ASC) and (3) no ASC (organic system without agroecological services crops). The qualitative response of the DEXi-met model suggested that the treatments with ASC returned similar total energy outputs, indicating the positive effect of this agroecological practice. Therefore, the authors pointed out that the ECO system could contribute to building up more complex agroecosystems, increasing both the resilience and biodiversity in organic agriculture.

Sofo and Ricciuti [21] presented a standardized method, named Biolog® EcoPlates™ (Biolog Inc., Hayward, CA, USA), to analyze the functional diversity of bacterial communities by means of measuring their ability to oxidize carbon substrates. Specifically, they reported a detailed methodological procedure which was easy to follow and reproduce (i.e., a detailed and step-by-step description of soil preparation, dispersion and dilution,

optimal bacterial density of the inoculum and so on), based on previous data and experiments. To this aim, a case study of the soils of a Mediterranean olive orchard was reported for comparing the methods and justifying the methodological protocol proposed. The authors emphasized that the results of this methodological paper could be important for correctly evaluating and comparing the microbiological fertility of the soils managed with sustainable, conservation-based practices or conventional, non-conservation-based ones.

Topic II comprises four papers. Zhang et al. [17] evaluated the impact of differentlevels of crude oil on the soil properties, physiological and chemical parameters of maize leaves and the phenanthrene content in the leaf. The results showed that (1) the soil water content significantly increased as the total petroleum hydrocarbons increased, and the soil electrical conductivity significantly increased compared with the control, and (2) the stomatal length and density, leaf K and Na contents and phenanthrene leaf concentration decreased in the contaminated soil compared with the reference group. The authors concluded that the stomata structure of maize could be influenced by crude oil, thus possibly controlling the accumulation of polycyclic aromatic hydrocarbons in aerial tissues. Practical implications for contaminated soil management were also hypothesized. Indeed, because stomata can play a key role in the aerial uptake of pollutants, regulation of stomatal movement can be usefully applied in phytoremediation of contaminated soil.

Petruzzelli and Pedron [18] investigated the influence of the soil characteristics on the tungsten uptake by maize for three different soils—Histosol, Vertisol and Fluvisol—in the Mediterranean area. Tungsten is largely used in high-tech and military industries. Therefore, soils are increasingly enriched in this element, and the possible transfer in the food chain represents a current issue in agroenvironmental studies. The results showed that tungsten concentrations in the roots and shoots increased with the increasing soil pH and decreased with the increasing organic matter. Further investigations were suggested by the authors to fully understand the mechanisms of tungsten transfer from the soil to the plant and the corresponding environmental impact of this element on human health.

Cesco et al. [19] reviewed the issue of copper (Cu) accumulation in European vineyard soils due to massive application from different sources (among others, metal-contaminated sludges and Cu-based pesticides for crop defense against pathogens). Although Cu is an essential micronutrient for equilibrated crop growth, excessive Cu concentrations in soil may lead to severe symptoms of toxicity, which are beginning to be recognized in crops and, particularly, in acidic soils. After dealing with the state of the art on this topic, they listed and discussed several current challenges, including (1) copper effects on soil agrobiodiversity, (2) rhizosphere management (management of soil–root–microorganism interactions), (3) biotechnologies and breeding for a more resistant plant material and (4) smart viticulture.

Sayara and Sánchez [20] reviewed the topic of bioremediation of contaminated soils due to the accumulation of polycyclic aromatic hydrocarbons (PAHs). Specifically, they focused on the application of compost as a type of biostimulation or bioaugmentation treatment to facilitate the bioremediation of soils contaminated with PAHs. Moreover, the authors reviewed the use of composting as an ex situ bioremediation strategy for PAH-contaminated soils. After dealing with the properties and sources of PAHs, they discussed the bioremediation of PAH-contaminated soils by composting and bioaugmentation (reviewing the effects of PAH characteristics and concentrations, temperatures of treatment, organic co-substrate stability and the mixing ratio), also describing the consequences of PAH bioavailability and biodegradation pathways. The authors concluded that, although the use of compost and composting in several remediation strategies significantly improved the removal of PAHs in contaminated soils, further investigations are needed in this research field.

#### **3. Conclusions**

The 12 manuscripts presented in this Special Issue can contribute to improving our understanding of sustainable agriculture and soil conservation through the development and

adoption of specific agronomic practices able to preserve and enhance the physical, chemical and biological properties of soils. The derived improvement of agroecosystem functions goes in the direction desired by both scientists and farmers. The investigations presented considered various crops, three different countries and multiple pedoclimatic conditions.

For the reader's convenience, the collected contributions were summarized in two main groups to provide results on the effects of soil use and management (topic I) and on the impact of some soil contaminants (e.g., crude oil, tungsten, copper and polycyclic aromatic hydrocarbons) (topic II).

The common conclusions drawn by the authors of this SI agree on the importance of both improving the available options for soil conservation and developing, implementing and sharing new techniques for the sustainable use of the soil.

**Funding:** This research received no external funding.

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

#### **References**


## *Review* **A Smart and Sustainable Future for Viticulture Is Rooted in Soil: How to Face Cu Toxicity**

**Stefano Cesco <sup>1</sup> , Youry Pii <sup>1</sup> , Luigimaria Borruso <sup>1</sup> , Guido Orzes 1,2 , Paolo Lugli <sup>1</sup> , Fabrizio Mazzetto 1,2 , Giulio Genova 1,3 , Marco Signorini <sup>1</sup> , Gustavo Brunetto <sup>4</sup> , Roberto Terzano <sup>5</sup> , Gianpiero Vigani <sup>6</sup> and Tanja Mimmo 1,2,\***


**Abstract:** In recent decades, agriculture has faced the fundamental challenge of needing to increase food production and quality in order to meet the requirements of a growing global population. Similarly, viticulture has also been undergoing change. Several countries are reducing their vineyard areas, and several others are increasing them. In addition, viticulture is moving towards higher altitudes and latitudes due to climate change. Furthermore, global warming is also exacerbating the incidence of fungal diseases in vineyards, forcing farmers to apply agrochemicals to preserve production yields and quality. The repeated application of copper (Cu)-based fungicides in conventional and organic farming has caused a stepwise accumulation of Cu in vineyard soils, posing environmental and toxicological threats. High Cu concentrations in soils can have multiple impacts on agricultural systems. In fact, it can (i) alter the chemical-physical properties of soils, thus compromising their fertility; (ii) induce toxicity phenomena in plants, producing detrimental effects on growth and productivity; and (iii) affect the microbial biodiversity of soils, thereby influencing some microbial-driven soil processes. However, several indirect (e.g., management of rhizosphere processes through intercropping and/or fertilization strategies) and direct (e.g., exploitation of vine resistant genotypes) strategies have been proposed to restrain Cu accumulation in soils. Furthermore, the application of precision and smart viticulture paradigms and their related technologies could allow a timely, localized and balanced distribution of agrochemicals to achieve the required goals. The present review highlights the necessity of applying multidisciplinary approaches to meet the requisites of sustainability demanded of modern viticulture.

**Keywords:** copper; rhizosphere; smart agriculture; microbes; vineyard

#### **1. Introduction**

In recent decades, the food demand has significantly increased in terms of quality and quantity due to the increase in the global population, which has now reached almost 8 billion people [1]. The growth of the per capita income in developing countries has undoubtedly played a decisive role [2,3]. On the other hand, the arable soil surface has decreased due to soil degradation and the impact of climate change [4]. Indeed,

**Citation:** Cesco, S.; Pii, Y.; Borruso, L.; Orzes, G.; Lugli, P.; Mazzetto, F.; Genova, G.; Signorini, M.; Brunetto, G.; Terzano, R.; et al. A Smart and Sustainable Future for Viticulture Is Rooted in Soil: How to Face Cu Toxicity. *Appl. Sci.* **2021**, *11*, 907. https://doi.org/10.3390/app11030907

Received: 18 November 2020 Accepted: 13 January 2021 Published: 20 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

we are experiencing a global annual mean warming of 1 ◦C above the level of the preindustrialization period [5]. Moreover, the global non-renewable natural resources used for the production of fertilizers (e.g., rock phosphate) are limited and are expected to be consumed shortly [6]. All these trends highlight the limits and the vulnerability of the current model of economic and social growth based on mass production and consumption, including with respect to soil [7].

For this reason, the United Nations defined 17 Sustainable Development Goals (SDGs) and 169 target actions for the year 2030, to which we are all invited to proactively contribute. These include targets of no poverty, zero hunger, good health and well-being, life below water, life on land, and climate action. In this context, agriculture is called upon to contribute by providing primary production and to reduce its anthropogenic impacts on the environment.

A specific example is represented by copper (Cu) accumulation in agricultural soils due to the massive application of metal-contaminated sludges and/or Cu-based pesticides for crop defense plans against pathogens. While Cu is an essential micronutrient for equilibrated crop growth, excessive Cu concentration, due to its accumulation in soil, leads to severe symptoms of toxicity, which are beginning to be recognized in crops, particularly in acid soils [8,9]. The case of vineyard soils is the most emblematic. Therefore, specific strategies must be adopted to limit and/or mitigate this problem, particularly in the viticulture sector.

This review aims at shedding light on the soil vineyard Cu accumulation and its agricultural and environmental consequences, specifically at the European level. Some approaches to limiting Cu toxicity and the related challenges, considering recent technological innovations (rhizosphere management, biotechnologies and precision/smart viticulture), are discussed. A short description of the role of Cu as an essential nutrient for plants is also included. The review further highlights research gaps that urgently require further study and innovation in order to guarantee sustainable vineyard management, providing soil conditions that will enable quality viticulture in the long term.

#### *1.1. Cu in Agricultural Soils and Crops*

Copper is a trace element in soil–plant systems. The average Cu concentration in the Earth's crust is 60 mg kg−<sup>1</sup> and in soil, it typically ranges between 2 and 50 mg kg−<sup>1</sup> . Natural soils with a high content of clay minerals (e.g., *Vertisols*) or organic matter (e.g., *Spodosols*, *Histosols*) are usually characterized by a higher Cu content (up to 180 mg kg−<sup>1</sup> ) [10]. Anthropogenic atmospheric depositions (traffic and industry related) and agricultural materials (inorganic fertilizers, agrochemicals, sewage sludge, livestock manure, irrigation water, compost, etc.) can dramatically increase the Cu concentration in soil. Indeed, Cu concentrations higher than 1000 mg kg−<sup>1</sup> have been recorded in polluted soils, including in agricultural soils [11]. In particular, concentrations exceeding 1000 mg kg−<sup>1</sup> have been reported for vineyard soils in France and Brazil [9] (see also the following chapter).

The main form of Cu in the soil is Cu2+, typically bound to inorganic and organic ligands, forming both soluble and insoluble compounds. In the solid phase, Cu has a very high affinity for organic matter and manganese oxides, as well as for iron oxides and clay minerals. Because of its high affinity for soil colloids and especially organic matter, Cu shows low mobility in near-neutral soils and is mainly concentrated in the upper soil layers. For instance, in rhizosphere soil (i.e., at the root–soil–microorganism interface), Cu is complexed by low-molecular-weight organic compounds exuded by plants and microbes like carboxylic and phenolic acids, which play an essential role in both external and internal plant tolerance mechanisms [9]. In arable soils, Cu concentration in the soil solution lies in the range 1–300 µg L−<sup>1</sup> [12].

Concerning the plant acquisition process, roots can use the free ionic forms (Cu+/Cu2+), although the direct use of Cu complexes with different organic molecules present in the rhizosphere cannot be excluded. Indeed, the identification of an uptake system of Fephytosiderophore (i.e., root exudates with a high affinity for Fe) complexes also in di-

cots [13] further corroborates this hypothesis. Moreover, a reduction step (Cu2+ → Cu<sup>+</sup> ) is also very likely to occur before the root uptake of the nutrient [9]. Indeed, previous works demonstrated that Fe chelate reductase enzymes and their related gene family (*FRO*) responded to both Cu deficiency and toxicity in *Arabidopsis thaliana* and cucumber plants, respectively [14,15]. The metal is taken up through a specific transmembrane Cu Transporter Protein 1 (COPT1) [16]. However, the involvement of other transporters like Zn/Fe Permeases (ZIPs), Natural Resistance Associated Macrophage Proteins (NRAMP) and Iron Related Transporters (IRT) has also been speculated [17–20]. In this context, the extent of the Cu bioavailable fraction (i.e., the two free ionic Cu forms and the soluble metal complexes with organic ligands) strongly depends on soil properties, especially cation exchange capacity (CEC), pH values, organic matter and clay content [12].

Once taken up by the roots, Cu is loaded into the xylem vessels to be translocated towards the aboveground organs. Despite only a few pieces of information available, the Heavy Metal-transporting P-type ATPase 5 (OsHMA5) has been shown to be involved in the Cu xylem-loading in rice plants [21]. In this respect, the Cu concentration in the xylem sap and in shoots is remarkably limited when the expression of this gene is down-regulated. Moreover, in *A. thaliana*, Yellow Stripe Like 2 (YSL2) transporter was also found to play a pivotal role in the xylem loading of Cu complexed with nicotianamine [22]. Once loaded in the xylem, Cu is bound to an organic ligand, generally amino acids (like histidine and proline in *Brassica carinata*) or nicotianamine (e.g., in tomato and chicory) [23–25]. Once at the leaf level, before being taken up by the cells, the Cu2+ complexes need to be reduced to Cu<sup>+</sup> by FRO4 [26] and, only after its transmembrane transport, the metal is delivered to the different organelles (e.g., mitochondria and chloroplasts) through specific transporters [27].

In plant cells, Cu is involved in a plethora of biochemical and physiological processes, including photosynthesis, respiration, oxidative stress responses, cell wall metabolism, and hormone signaling [28,29]. In fact, Cu is an important cofactor in many proteins, such as plastocyanin, cytochrome c-oxidase and amino oxidase [28]. Generally, the Cu concentrations in plant tissues growing in uncontaminated soils range from 5 to 30 mg kg−<sup>1</sup> dry weight, depending on the type of plant, the growth stage, and the soil characteristics. Deficiency symptoms can appear at Cu concentrations lower than 5 mg kg−<sup>1</sup> dry weight, while leaf concentration higher than 20 mg kg−<sup>1</sup> dry weight may result in toxic effects affecting the whole plant development [30–32].

#### *1.2. Current Situation in Vineyards (Including Toxicity in Plants and Effects)*

Several studies have reported an accumulation of metals in agricultural soils due to common agricultural practices, such as manure fertilization and spraying for pest control [33–35]. For instance, Cu in apple orchards can have an average annual total increase ranging from 2.5 to 9 mg kg−<sup>1</sup> year−<sup>1</sup> deriving from the application of Cu-containing agrochemicals [36]. The accumulation of Cu through time can be influenced not only by the age, but also by other factors, such as soil organic matter content and pH [33]. As shown in Figure 1, vineyard soils are prone to Cu accumulation, since the practice of spraying Cu salts as fungicides dates back at least as far as 1761 [37]. After the introduction in 1885 of the mixture of Cu sulfate and lime, called "Bordeaux mixture" [38], its use became generalized on wine-producing farms. Because of this, the concentration of Cu in vineyards can be many times higher than the background values of forest soils in the same region, reaching up to 3000 mg kg−<sup>1</sup> [39,40]. Next, high Cu concentration differences can also be found inter-row and intra-row within the same vineyard [41].

‐

−

‐

‐ ‐ ‐

‐

‐

‐

−

−

‐

‐ ‐

‐ ‐ ‐

**Figure 1.** Cu concentration distribution (**a**) compared to the distribution of vineyards (**b**) in Europe [42]. Several high Cu concentration areas perfectly match viticulture areas.

‐ − ‐ ‐ ‐ ‐ Although much research is currently taking place to find alternatives to Cu, this element is still crucial to fighting fungal diseases in vines, especially in organic farming, limiting the supply of pure Cu to 28 kg ha−<sup>1</sup> over a period of 7 years [43–45]. For these reasons, the global trends of vine-cultivated areas and the use of Cu should be monitored, particularly in newly planted vineyards. Indeed, the shallow root system of newly planted vine plants might be more affected by high Cu levels than older plants, since Cu mainly accumulates in the upper soil layers (where the concentration of organic matter is higher). In fact, the first effects of Cu phytotoxicity are evident at the root level, with a clear decrease of root elongation, abnormal root branching, thickening and dark coloration [46]. As a consequence of this altered root development, the roots' ability to take up water and other nutrients is clearly impaired [32,47]. In this regard, it is worth mention the negative effect of Cu toxicity on P acquisition by roots [48,49], highlighting the interactions occurring between the different nutrient acquisition mechanisms. Critical Cu concentrations in roots range from 100 to 400 mg kg−<sup>1</sup> dry weight [10,46]. Like other heavy metals, an excess of Cu may also cause the generation of toxic reactive oxygen species (ROS), which, in turn, can damage several important biomolecules like DNA, proteins and lipids [50]. Usually,

−

‐

‐

plants exhibit an excellent translocation barrier of excessive Cu between roots and shoots. This limits the risk of Cu accumulation in edible plant tissues that exceeds toxic levels to livestock and humans. However, toxic concentration levels, as well as symptoms of Cu toxicity, have also been well described for the aerial part of plants [30–32]. Typical cases of Cu toxicity at the root and shoot levels of vine plants (*Vitis vinifera* L., *cv Glera*) are presented in Figure 2.

**Figure 2.** Symptoms of Cu toxicity in vine plants grown in a soil with 814 mg kg−<sup>1</sup> Cu (DTPA extractable fraction: 280 mg kg−<sup>1</sup> Cu): (**D**–**F**) canopy, (**C**) leaves and (**H**) roots of grafted plants of *Vitis vinifera* L. (*cv. Glera*); (**G**) canopy of the rootstock SO4 (*V. berlandieri x V. riparia*). As control, canopy (**A**) and leaf (**B**) of healthy vine plants (*Vitis vinifera* L. *cv Glera*).

‐ ‐ ‐ ‐ Considering the soil surface dedicated to viticulture, in 2019, worldwide vine cultivation covered 7.4 million hectares [51], and it has been considered stable in recent years, with a global balance given by heterogeneous evolution in different regions of the world [52]. In Europe (Figure 1), the vine-cultivated area is stable, being balanced between the European Union "grubbing up" program (i.e., the replacement of vine plants with other agricultural crops, [53]) on the one hand, and the possibility of Member States authorizing the planting of up to 1% of the vineyard surface already planted each year [54], on the other hand. Globally, several countries are increasing their vineyard area, namely Russia, New Zealand, Peru, and China, to name a few. Meanwhile, countries like USA, Argentina,

‐

‐

Chile are showing a decline in vineyard area [51]. In addition to wine regulations and the global market [55], Climate Change is playing a fundamental role in directing the spatial distribution of vines, as increasing temperatures are leading to vineyards higher in elevation and at higher latitudes [56,57]. Climate change also affects the incidence of fungal infections in vineyards, requiring increasing use of Cu-based treatments [9]. These changes will then affect the chemical and biological quality of soils of newly planted vines with increasingly higher Cu concentrations. In this regard, it is interesting to note that the annual supply of Cu for pathogen defense plans can range from 1–2 kg ha−<sup>1</sup> in Europe to 30 kg ha−<sup>1</sup> in southern Brazil [58] as a function of the infectious pressure of the pathogens. Furthermore, new crops growing where there used to be vineyards might be affected by high Cu concentrations and exhibit symptoms of toxicity.

European countries have different approaches to defining threshold levels, as there is large heterogeneity in the legal systems, the chemical analysis used, the target organisms for toxicity, and background-values and how these are defined. In addition, other soil properties (mainly pH and SOM) and the precipitations/humidity of wine areas can be used as additional data to determine the final threshold level since the need for applying higher quantities of active principle can influence the decision-taking process [59–61]. In Italy, for instance, the management of polluted sites is regulated by the law D. Lgs. 152/2006, and a recent Decree established 200 mg kg−<sup>1</sup> of Cu as the contamination limit for agricultural soils [62]. A literature review showed a high degree of heterogeneity and a lack of standard approaches and thresholds for heavy metal risk assessment [59]. For example, the Austrian standard ON S 2088-2 defines a two-step evaluation in which, firstly, soil metal content is measured (NH4NO3-extract). If certain trigger values are reached, then the possible bioavailability is evaluated in consideration of the metal content and other soil parameters, such as pH and soil organic matter, thus resulting in a final contamination evaluation [59]. The Czech Republic provides indication limits on plant growth inhibition [63] in consideration of the metal content extracted with aqua regia (HCl and HNO<sup>3</sup> at a ratio of 3:1) or with 1 mol L−<sup>1</sup> of NH4NO3. Both extractions must be performed if the limit values for the specific metal are defined. At the global level, the lack of a standardized approach on Cu content thresholds in soils should be addressed in the future both by the scientific community and in public policy. Several soil extraction methods should be evaluated, considering both total and available Cu contents in order to reach a robust decision scheme.

#### **2. Current Challenges**

As described in the previous paragraphs, a deeper understanding of the Cu accumulation phenomenon in vineyard soil with an all-encompassing approach, including the main environmental and agricultural risks connected, appears necessary. Concurrently, the identification of approaches and the setup of cultivation practices aimed at mitigating the problem are urgent. Therefore, the following paragraphs describe some examples considering recent technological innovations like rhizosphere management, biotechnologies applied to vine plant breeding and smart viticulture.

#### *2.1. Copper Effects on Soil Agrobiodiversity*

Soil contains a vast underground world inhabited by many fungal and bacterial species. It is estimated that 1 g of soil contains a number of bacteria ranging between 100,000 to 1,000,000, and, within this fraction, there may exist different micro-niches in which several ecotypes may live [64,65]. Although only a small fraction of the fungal and bacterial diversity is known, certain fungal and bacterial groups play a crucial role in several agroecosystem services, including plant-growth promotion and cycling of nutrients (e.g., nitrogen, N, and phosphorus, P) [66–68]. However, the accumulation and the excessive availability of heavy metals such as Cu can impact the soil microbial biomass and activity, thus affecting some microbial-driven soil processes, including N fixation and mineralization [69]. Moreover, Cu has been demonstrated to impact microbial

communities' composition and functionality, with higher effects on bacterial than fungal communities [70–74]. Fungi and bacteria have evolved different strategies against high Cu concentrations in soils (Figure 3). For instance, some filamentous fungal and yeasts species act against the overload soil Cu, directly at the cytosol level [75]. Herein, the sequestration of the free Cu ions occurs due to metallothioneins, a family of cysteine-rich proteins (e.g., Cup1 and Crs5) that chelate Cu [76–78]. Furthermore, when the metal in the cytosol is in excess, fungal cells act through organellar compartmentation of the metal, mainly directing it to the vacuole. Concerning this latter strategy, no specific Cu transporter catalyzing the transmembrane transport into the vacuole has yet been identified [79]. In addition, cytoplasmic chaperones' activity is essential for mitigating Cu stress in the cytosol. Indeed, while the chaperone Atx1 delivers Cu to the trans-Golgi network, Ccs1 delivers Cu to the Cu-Zn Superoxide dismutase Sod1 [80,81].

‐ ‐ ‐ ‐ **Figure 3.** Exposure of soil microorganisms to Cu in vineyards. (**a**) Selection of Cu-resistant and antibiotic-resistant genes due to the Cu application (**b**) Schematic representation of the bacterial and fungal strategies for coping with excessive Cu inputs;

‐ ‐

‐

‐ ‐

‐

‐ ‐ ‐

‐

‐

‐

In contrast, in bacteria, the detoxification strategy relies mainly on a Cu export system based on the functionality of two different mechanisms for Gram-positive and Gramnegative bacteria [82]. In Gram-positive bacteria, proteins belonging to the Cop family are at the basis of cytosolic Cu disposition. The system includes a Cu chaperone (CopZ) that delivers Cu first to CopY (copper responsive repressor), activating the transcription of *copA* and *copB* (genes encoding for ATPase pumps), which extrude Cu<sup>+</sup> and Cu2+ from the cell, respectively. In Gram-negative bacteria, the presence of periplasmic space determines the evolution of specific mechanisms for extruding Cu in the extracellular space. Specifically, the Cop-protein system is located on the inner membrane and in the inner cell, and its activation is mediated by the MerR-type Cu-responsive transcriptional activator CueR [83]. Additionally, while Gram-positive bacteria mainly use the Cop system, Gram-negative bacteria exploit a trans-envelope extrusion system (CusABC ATPase pumps system) [82]. Additionally, the periplasmic multicopper oxidases (PcoA, CueO), oxidizing Cu<sup>+</sup> to Cu2+ and activating the production of siderophores to sequester the cytoplasmic Cu, play a crucial role in the Cu detoxification process [84,85].

Moreover, the Cu accumulation in agricultural soils other than bacterial Cu resistance can also promote the occurrence of antibiotic resistance, with these traits being linked together [86–90] (Figure 3). Therefore, at the molecular level, the Cu-resistant and antibioticresistant genes can be reasonably co-selected. Thus, they can commonly occur on the same mobile genetic element (i.e., plasmid, integron or transposon) [91,92]. Additionally, the bacteria cells can have a cross-resistance molecular mechanism in which the same gene confers resistance to different antimicrobial agents, such as antibiotics and Cu [91]. Indeed, the repeated application of Cu to agricultural soils over the years, given that the metal is not degradable, can dramatically enhance the frequency and dissemination of Cu-resistant and antibiotic-resistant genes [92]. In this context, Cu-contaminated vineyard soils could even be considered a reservoir of antibiotic-resistance traits that can be transferred among the different bacterial species via vertical and horizontal gene transfer (e.g., plasmid-mediated conjugation, integrons) [36]. However, more detailed studies should be carried out to elucidate the possible long-term impact of Cu on the bacterial communities inhabiting agricultural soils linked to the antibiotic resistance. Indeed, there is a general lack of knowledge about the potential route for the transmission of antibiotic-resistant bacteria from soil to crops, animals and humans. Specifically, it will be crucial to investigate the potential health risk assessment to prevent further development of pathogenic-resistant bacteria in agricultural ecosystems.

#### *2.2. Soil Management Considering the Root–Microorganism Interactions (Rhizosphere Management)*

Considering the ever-increasing incidence of soils contaminated with Cu, it is evident that rhizosphere dynamics must be taken into account when setting up and applying soil management practices to mitigate the Cu stress in these agricultural soils. In fact, the biogeochemical cycles of the nutrients completely differ at the soil–root interface compared to non-rooted soil (i.e., bulk soil). Indeed, the highly dynamic interactions between roots, soil and microorganisms occurring in the rhizosphere control not only the speciation and the consequent availability of the nutrients (and, in turn, the extent of their root acquisition), but also determine the level of availability/toxicity of heavy metals like Cu. For instance, in the hotspot that is the rhizosphere, pH value modifications and the release of low- and high-molecular-weight organic compounds by both roots and microorganisms are the main drivers of Cu availability.

Concerning soil pH values, root-induced changes of this chemical parameter are mainly triggered by plant nutrient uptake. This process is often coupled with an active proton extrusion catalyzed by the plasma membrane H<sup>+</sup> -ATPase [28]. Furthermore, an unbalanced cation/anion uptake of nutrients might lead to either acidification (i.e., due to an excessive uptake of cations) or alkalization (i.e., due to an excessive uptake of anions) of the rhizosphere. These imbalances might occur naturally due to different plant requirements, depending on the physiological status of the plants, microbial interactions

and/or competition with neighboring plants, as well as potentially being caused by the unbalanced availability of some nutrients in the growth medium [93,94]. Indeed, in vineyards, when cover crops are co-cultivated with other plants like vine plants, nutrient competition might induce an unbalanced uptake of elements with consequent root-induced changes of the pH. As is to be expected, this phenomenon impacts Cu mobility/availability and, consequently, its toxicity. On the other hand, excessive concentrations of metals like Cu might also be themselves responsible for an unbalanced nutrient uptake. Recent studies have highlighted that Cu stress can induce synergism and/or antagonisms with many other nutrients like P, Fe, Zn and Mn, causing further root-induced pH changes both in annual grass species and in perennial plants like grapevine [15,95]. The complexity is further exacerbated by the fact that rhizosphere effects (like acidification and alkalization processes) are not always constant over time [96] and are strongly plant species and genotype/rootstock dependent [15,95,96].

As mentioned earlier, besides pH values, rhizosphere organic compounds are the main drivers shaping Cu speciation in agricultural soils. Indeed, specific compound classes such as phenolic compounds and carboxylic acids play a fundamental role in both internal and external Cu tolerance strategies of plants [97]. Recent studies have shown that perennial plants like grapevines also trigger their root exudation when exposed to excessive Cu concentrations, yet the phenomenon strongly depends on the type of the rootstock and the growing conditions [15]. Indeed, when grapevines are intercropped with annual grass species like oat, the exudation pattern is completely modified. Furthermore, the enhanced exudation reduces the Cu accumulation in grapevine plants, making this a valuable agronomic strategy for mitigating Cu stress [15]. It is important to highlight that the alleviating effect is highly plant species dependent, and thus the intercropping approach should be evaluated for each specific rootstock. In addition, the intermingling of roots in intercropping systems might induce a competition between plant species leading on one side to a modified exudation profile and on the other side to an unbalanced nutrient uptake ultimately affecting the effect of Cu stress.

In addition, soil management also comprises fertilization strategies, which also strongly affect rhizosphere dynamics. For instance, the source of N supply (NO<sup>3</sup> <sup>−</sup> or NH<sup>4</sup> + ) affects the pH at the soil–root interface both in annual and perennial plant species [28]. Nitrate-based fertilizers induce an enhanced anion root uptake with proton consumption and a consequent alkalization of the rhizosphere. On the other hand, ammonium-based fertilizers trigger enhanced proton release by roots, resulting in significant rhizosphere acidification. Even though nitrate usually predominates over ammonium in agricultural soils due to the very rapid microbial nitrification processes, preferential ammonium uptake might still occur in acid soils, at low soil temperatures, or shortly after the application of ammonium fertilizers, organic fertilizers and/or nitrification inhibitors [28]. Soil management should thus comprise fertilization-induced root activities, since Cu availability in soil is strongly pH dependent and might undergo mobilization processes, exacerbating its toxicity. Moreover, an adequate nutritional balance for a crop seems to be decisive in mitigating the effects of Cu toxicity when the metal is already present in relevant quantities. Vineyards in calcareous soils represent a clear example. In fact, in this case, the edaphic conditions (essentially soil pH values) should theoretically restrain the availability of Cu to the plants. However, if the Fe availability for crops becomes limited as a consequence of the chemical–mineralogical nature of the soil, the activation of the plant's adaptive responses to this deficiency (including the acidification of the rhizosphere [98]) could transform the Cu problem from potential to real.

The complexity of rhizosphere dynamics highlights that soil management, and particularly the management of nutrients, in agricultural agroecosystems such as vineyards needs a comprehensive overview and an appropriate knowledge of the soil–plant system and its dynamics. Indeed, such knowledge could be crucial in counteracting the negative effects of toxic concentrations of heavy metals like Cu and/or exploiting beneficial synergism between elements to benefit the agricultural production. Future studies should therefore

focus on a holistic approach that also includes the third main actor involved in rhizosphere processes, i.e., microorganisms. Even though microorganisms are known to increase plant growth and stress resilience, little is known about the effect of Cu stress on their composition and functionality in the rhizosphere in real agroecosystems such as vineyards. To date, most studies have been lab-based and recreate artificial metal stress conditions; thus, they still provide an insufficient understanding of the microbiome functioning and the mechanisms of plant–microbiome–soil interactions. Large-scale experiments at the field level are thus essential in order to give a complete overview of the effect of Cu on the complex interplay between soil, plants and microorganisms involved in the biogeochemical cycles of nutrients at the agroecosystem level.

#### *2.3. Biotechnologies and Breeding for a More Resistant Plant Material*

Grapevine cultivation and wine production worldwide face challenges posed by the high pressure of fungal and fungal-like diseases, causing production losses [99]. The most common and severe problems in viticulture include those presented by downy mildew and powdery mildew infections, caused by *Plasmopara viticola* and *Erisyphe necator*, respectively [100]. Both pathogens are commonly controlled by the application of fungicides in vineyards, mainly based on Cu [101]. However, the implementation of innovative strategies, such as the use of pathogen-resistant vine genotypes, is highly desirable with the aim of achieving a more sustainable viticulture. For this reason, breeding programs have been implemented in order to transfer the so-called resistance (*R*) genes from resistant *Vitis* species to sensitive *V. vinifera* varieties [102].

The production of resistant grapevine hybrids began in the 19th Century, yet the first varieties, obtained from traditional breeding carried a high percentage of non-*V. vinifera* genetic material and presented undesirable "foxy" flavors in the resulting wine, a distinctive feature of *V. labrusca* [101]. In the post-genomic era, new techniques, such as marker-assisted backcrossing, marker-assisted background selection and marker-assisted pyramidization, contributed to the development of breeding strategies allowing the precise introgression of wild alleles in susceptible genomes, while reducing the undesired residual genetic material [103]. At present, up to 27 Quantitative Trait Loci (QTL), known as *Rpv* genes, are associated with the resistance to *P. viticola*, while up to 13 QTL (*Ren* and *Run* genes) are related to resistance against *E. necator* [102]. Following infections, *Rpv* genes can induce different defense responses (e.g., hypersensitive response, the accumulation of phenolic compounds in the infected tissues, the accumulation of callose and lignin, the synthesis of phytoalexins, the induction of either cell necrosis or peroxidase activity; [104] and the references therein), while *Ren* and *Run* loci have been shown to encode small gene families of receptor-like proteins that are thought to directly or indirectly interact with pathogen-specific effectors. This interaction elicits a signaling cascade that leads to the transcriptional reprogramming of the host plant and the expression of plant defense genes [105]. Despite the potential economic and environmental benefits achievable, the diffusion of resistant hybrids is being prevented by the wine-making industry, which preferentially chooses traditional genotypes, despite their being susceptible, over resistant ones, mainly because of current regulations [102,106]. A remarkable opportunity to overcome the issues related to traditional breeding might be offered by the genome editing approach, directed towards the susceptibility (*S*) genes encoded in the *V. vinifera* genome [107]. The functionality of *S* genes is required for successful infections by both *E. necator* [108,109] and *P. viticola* [106,110]. Indeed, the targeted inactivation of *S* genes in *V. vinifera* plants might potentially lead to the generation of resistant clones, while still preserving the genetic identity of the parental genotypes, which would be highly appreciated by the wine producing industry [106].

Interestingly, these new varieties can be defined as being tolerant to the pathogens; however, in cases of very high infection pressure, the endogenous plant responses might not be sufficient for counteracting disease, making the application of agrochemicals like Cu necessary for controlling and/or contrasting pathogen diffusion within the canopy

‐

(Figure 4). Indeed, despite not completely avoiding the application of fungicides, resistant varieties might contribute to strongly decreasing the use of agrochemicals in open fields.

‐ ‐

‐

‐

‐

‐

‐

‐ ‐ ‐ ‐ **Figure 4.** Traditional susceptible varieties vs. resistant varieties in viticulture. The protection of traditional susceptible varieties from pathogen infections requires the repeated application of Cu-based fungicides, which, through drift phenomena and rain wash-offs, can cause a consistent accumulation of Cu in vineyard soils. The increased Cu concentration in soil can, on one hand, induce vine plants to modify their exudation to cope with Cu toxicity, and, on the other hand, promote an alteration of soil microbial biodiversity. The combination of modified exudation patterns with an altered rhizosphere microbial community might impact the biogeochemical cycles of mineral nutrients, thus modifying their availability to plants. In resistant varieties, the protection against pathogens requires lower application rates of agrochemicals, thus also reducing the input of Cu. Nonetheless, the expression of the resistance gene(s) might alter the physiology of vine plants at transcriptomic, proteomic and metabolomic levels, in this case also causing a modification in the exudation pattern. In addition, evidence shows that resistant varieties might require specific fertilization strategies in order to fully present the resistance traits. Again, an altered exudation pattern combined with the input, albeit minimal of Cu and fertilizers may impact the rhizosphere dynamics, affecting both soil microbial biodiversity and the biogeochemical cycles of mineral elements.

Nevertheless, both the introgression of foreign genetic material (such as for traditional breeding) and the silencing of genes (such as for the genome editing approach) could potentially alter the physiology of plants. The modification of vine plants' genetic material might produce, for instance, organisms that could require particular edaphic conditions (e.g., specific fertilization practices and/or specific provision of micro/macronutrients) in order to express disease resistance at optimum levels (Figure 4). In fact, specific nutrient fluxes within the leaf blade seem to be the basis for the expression of the pathogen response in these resistant varieties [111]. In addition, transcriptional reprogramming, as in the case of either *Rpv* or *Ren/Run* loci, can alter the metabolome [112,113], which might mirror modifications in qualitative and quantitative root exudation patterns. Root rhizodeposition plays a crucial role in the transformations and fluxes of nutrients from soil to plant [114], considerably affecting the extent of their acquisition by plants, and therefore, the crop yield. Moreover, root exudates represent, at least in part, the mechanism by which plants can recruit soil microbiota and shape the rhizosphere microbial community [115]. As mentioned earlier, this could either directly or indirectly impact the structure, the diversity and the functionality of soil and rhizosphere microbial communities, which play a pivotal

role in the processes occurring below ground (Figure 3) [68,115–117]. In this regard, the accidental selection of pathogenic microbes, which could represent a future challenge for viticulture in terms of new plant diseases, cannot be excluded. Therefore, the suitability assessment of these new plant materials should not be limited to those traits that are strictly connected with the expression of the pathogen resistance genes for the purpose of overcoming old (i.e., the classical and well-recognized vine diseases) challenges. Indeed, it is necessary to evaluate resistant varieties also with respect to putative new challenges (i.e., new pathogens and the related diseases) in order to establish the sustainability of these viticultural practices in both the short and the long term. Furthermore, the availability of rootstocks with traits of tolerance to Cu toxicity could be advantageous in viticultural areas. Finally, the increase in phenotype and genotype data of the scions and the rootstocks could significantly benefit from the application of machine learning to accelerate the crop breeding process [118] in press).

#### *2.4. Smart Viticulture*

Cu accumulation in soil can be restrained through direct (e.g., using more resistant plant species, as described in the previous paragraph) and indirect actions (e.g., decreasing the availability in the rhizosphere and/or limiting the quantity of Cu that reaches the soil via new applications). Since the resistant varieties of vine plants currently available are essentially tolerant, and therefore a certain level of plant defense using pesticides is essential, indirect actions are still important not only for the traditional varieties, but also for the new ones. In this respect, the main goal is to avoid drift phenomena during the application of pesticides and to limit Cu accumulation in the soil. The Cu applied to the vine canopy will indeed reach the soil due to rain, foliar irrigation and/or autumn leaf fall.

In the last two decades, a wide range of technologies have been increasingly applied in viticulture, such as monitoring technologies, decision support systems and operating technologies. **Monitoring technologies** include geolocation, remote sensing (using satellites, aircrafts, and unmanned aerial vehicle/drone images), and proximal sensing during production (using wireless sensor networks to measure soil moisture, leaf wetness, grape temperature, sap flow, etc.) and harvesting phases (using volumetric grape sensors and optical sensors for grape "quality"). **Decision support systems** (DSS) make it possible to consider the spatial variability (i.e., the variability of soil properties, landscape features, crop stresses, yield and quality in the different areas of a field [119]) in process optimization (irrigation, fertilization, or chemical treatments) and harvesting. Finally, operating technologies include **variable-rate technologies (VRT)** and **agricultural robots**, which make it possible to implement the activities defined with the help of the DSS, while also limiting drift phenomena [120].

These technologies were proposed in the 1990s, initially for their ability to manage crops according to site-specific logic, which makes it possible to deal with the spatial variability of soil, nutrient and phytosanitary conditions. They also include geo-referencing, which allows producers to micro-manage soil and plant processes within small areas of a single field. Despite positive expectations, these techniques' diffusion is still relatively limited, as several authors have reported [121–125]. In several European countries, the use of ITC in agriculture remains at less than 10% of farms, and, of these, only 50% of farms that complete site-specific applications. The current circumstances are considered to be very favorable for the growth of the sector, not only for the new generation, with habits more rooted in the use of digital technologies, but also for the "cultural dragging" caused by the so-called "Industry 4.0" revolution [126–128]. In fact, this revolution has brought significant attention to the new technological frontiers connected to the "Internet of Things" (IoT), Big Data, Cloud Computing, hyper-connectivity, cybernetics and augmented reality. As a consequence of this technological trend and its effects, the terms *agriculture 4.0* and *smart agriculture* are coming to be preferred to the term *precision agriculture* [129].

To the best of our knowledge, there are no studies measuring the exact impact of precision or smart viticulture techniques on the limitation of Cu accumulation/toxicity

at the soil level. However, various papers have highlighted that the goal of these techniques is to improve not only the yield and quality of grape production, but also its environmental sustainability (i.e., reducing the chemical treatments) [120,130]. Indeed, these techniques have shown that it is possible to reduce the use of Cu-based fungicides by 25–50% [131–133]. Therefore, it can be expected that the contribution to the limitation of Cu accumulation/toxicity at the soil level will be remarkable. Thus, a set of examples of *precision* and *smart viticulture* technologies (which might be used by practitioners and scholars as an initial reference point for implementation projects and precise impact analyses) is discussed in the following paragraphs. Other examples can be found in the detailed reviews of literature and field applications of precision viticulture [120] or precision agriculture [134–136].

While soil protection has mostly been studied using a mono-disciplinary approach (e.g., soil chemistry, plant defense), the use of precision/smart viticulture technologies in combination with an interdisciplinary approach involving technology experts (from the fields of sensors, electronics, computer science, and agricultural machinery), domainexperts (e.g., from the fields of soil chemistry, plant pathology, genetics, and agronomy) and sustainability experts (or Lifecycle assessment experts) is increasingly needed.

In particular, *precision* and *smart viticulture* technologies could contribute positively to the following aspects: (1) crop monitoring, i.e., using optical sensors on vehicles (e.g., tractors, autonomous vehicles, or drones) to perform a sort of "early digital scouting" to be carried out promptly to keep the phytosanitary status of wide-crop areas under control [137,138]; (2) operational monitoring, i.e., using the so-called FORK systems (Field Operation Register Keeping) [129,139,140] which allow the detection and automatic recording of how a field operation is carried out; (3) implementation of site-specific techniques with retrofit components, i.e., adaptable to existing farm machinery, thus avoiding the need to make large investments in new equipment (see Figure 5). The first two aspects concern information management actions, necessary for medium- and long-term quality decision-making processes at the farm, based on targeted information. The third aspect allows a direct intervention with immediate control effects, especially if the retrofit device is also equipped with sensors capable of locally estimating the volumes of canopy to be treated and then enabling the adjustment of the doses to be distributed accordingly.

Sensors installed both on agricultural machines and within the vineyards (both at the plant and soil level) allow the collection and monitoring of data (e.g., the color, shape and dimensions of leaves and grapes, temperature and weather conditions, soil pH, and the presence of pathogens [120,130,131,141]). These data can then be processed (cleaned, filtered, aggregated, represented and archived) and analyzed to extract the relevant information. These analyses might be carried out by a DSS at the farm or at external servers accessible over the Internet (i.e., in the Cloud) [133,141]. Moreover, they might be automated through Artificial Intelligence and/or neural networks algorithms, and they might also be used to simulate/prevent future behaviors/events (e.g., a digital twin of the vineyards can be created where the impact of different interventions/strategies can be simulated) [142]. The data analyses might also combine data from multiple farms/vineyards, as well as other big (information assets characterized by such a high volume, velocity and variety that it requires specific technology and analytical methods for its transformation into value [143]) data from several sources [144,145]. Finally, the results obtained can lead to an action, such as process optimization (modification of irrigation, fertilization, or chemical treatments), or to a report for internal or external users [129,131,133]. Such actions can also be implemented in a more precise and/or automatic way through straddle tractors/machines, autonomous vehicles/robots, and drones [120,146]. Finally, the application of these technologies might help the winegrower to spray, fertilize or irrigate (1) when it is needed; (2) exactly where it is needed (also limiting the drift); and (3) at the quantity needed in the different areas of the field and times (i.e., considering the horizontal, vertical, and spatial variability), thus potentially minimizing the Cu accumulation process in the soil. The technologies presented above are (mostly) already available on the market

and have been adopted by some vineyards. An interesting example of this is PlantCTTM (formerly known as Smartvineyard), a system that continuously monitors the plants and the environment using a wide range of sensors (i.e., spore and pest detectors, as well as leaf surface, light, humidity, and temperature, wind, solar radiation, precipitation and soil sensors) and suggests to the winegrower the interventions/actions that should be implemented (https://plantct.com/). Another example—more focused on the topic of this paper—is Coptimizer, a model-driven DSS that might help to optimize and track the use of Cu-based fungicides in viticulture. Other similar examples include VineSens [141], GRover [147], and FeelGridTM (https://www.feelgrid.com/). For a detailed review of smartphone applications for smart agriculture, see [148]. There are then several companies producing variable-rate (or smart) straddle tractors/machines and spraying drone systems (also) for vineyard applications (e.g., Pellenc, New Holland, Durand-Wayland, Tecnovit, Fly Dragon Drone Tech., Chouette). Some of these systems have been presented in detail by [120]. ‐ ‐ ‐ ‐ ‐ ‐ ‐

‐ ‐

‐

‐ ‐

‐

‐

‐ ‐ ‐ ‐ ‐ ‐

‐

‐ ‐ ‐ ‐ **Figure 5.** Design of a possible retrofit device to be applied to sprayers already present at the vineyard to enable them to operate according to site-specific logics. The retrofit application concerns a vertical spraying tower equipped with a set of rotating supports for independent nozzle selection and activation. The control is supported by a set of volumetric ultrasonic sensors able to detect the size of the canopy to be processed. The solution even includes a Field Operation Register Keeping (FORK) system based on an identification device enabling the tractor to automatically detect the identity of the coupled machine through a simple Radio Frequency (RF) system (T-MO: implement transmitting code; A-RF: tractor receiving antenna; identification is enabled only on distances <20 m). On board data transmission is performed via Controller Area Network (BUS-CAN), supposing this is already available on-board. In case of older tractors (no BUS-CAN), a direct wire connection can be easily established.

The necessary investments and/or operating costs required by the abovementioned systems—as well as more generally for the implementation smart farming solutions—are, however, rather high, particularly for small vineyards [149,150]. Furthermore, specific know-how is often needed to properly use them [149]. The Industry 3.0 phase (indicatively occurring between 1970 and 2010) was characterized by digital culture and technologies in the management of industrial production processes. On the other hand, the Agriculture 3.0 phase (occurring between 1980 and 2000) introduced technological innovations mainly focused on the improvement of the quality of mechanization in the fields of electronics, ergonomics and safety, but with limited results in the field of information technology. The "digital gap" in agriculture compared with the industrial sector is particularly relevant. Therefore, the key challenge is to adopt these systems and technologies for economic sustainability (the retrofit idea mentioned above could, for instance, be an interesting solution) and to explain their potential benefits to winegrowers. Another possible idea could be

the creation of some consortia or cooperative forms to cover the high investments/costs. In any case, with the application of these strategies, it is possible to significantly limit the supply of Cu to the vineyard, thus limiting the metal accumulation process in the wine-growing context.

However, future research is needed on the smart viticulture topic at different levels to contribute to the general challenge. First, while well-developed and robust solutions for specific purposes exist (e.g., grape yield monitors, soil and leaf monitors and precision sprayers), they still tend to work as separate silos. Future research should therefore try to develop visions as well as practical architectures/applications for an integrated smart (or digital) management of vineyards. Second, there is a strong need to identify the competencies needed by farmers (and by wine-growing consultants) in order to correctly and effectively use the technologies available and to consequently update the study programs. Finally, smart viticulture's impact on the three sustainability dimensions (i.e., environmental, social and economic) should be further investigated. [151] showed that both positive and negative impacts of Industry 4.0 technologies can be observed in the industrial field. This might also be applied in viticulture (and, more generally, in agriculture) where, to date, research has mostly focused on the environmental dimension [130,131].

#### **3. Conclusions**

The repeated application of copper (Cu)-based fungicides in conventional and organic farming has caused a significant increase in the total Cu contents in vineyard soils, posing agricultural and environmental threats. The goal of our review was to shed light on this issue and to discuss some approaches that might be adopted to limit Cu toxicity.

We presented the current Cu accumulation situation, with particular reference to the European context (see Figure 1), and discussed its potential risks in viticulture (see Figure 2), considering the chemical-physical-microbiological properties of soils, as well as the toxicity phenomena in plants. Furthermore, some approaches and setups of cultivation practices aimed at mitigating the Cu accumulation problem (namely, rhizosphere management, biotechnologies and breeding for more resistant plant material, as well as precision/smart viticulture techniques for a timely, localized and balanced distribution of agrochemicals) were discussed. We summarized the current literature for all these topics and highlighted a set of research gaps that urgently require future studies and innovation. In particular, research is needed to shed light on the possible long-term impacts of Cu on selecting antibiotic-resistant bacterial species of agricultural soils and on their potential threats to animal and human health. Future studies should then focus on vineyard rhizosphere dynamics more holistically, i.e., considering the microorganisms and the effect of Cu stress on their composition and functionality, as well as relying on large-scale experiments at the field level. Third, the assessment of "new" resistant plant material (scions) should not be limited to the traits strictly connected with the expression of the pathogen resistance genes for overcoming current vine diseases, but should also consider putative new pathogens and related diseases. From this perspective, particular attention should also be paid to rootstocks with tolerance traits for high Cu concentrations in soils. Finally, future research is needed (1) to develop visions as well as practical architectures/applications for an integrated smart (or digital) management of the vineyard; (2) to identify the competencies needed by farmers and wine-growing consultants to use the technologies available correctly and effectively; and (3) to shed empirical light on the impact of smart viticulture on the three sustainability pillars (i.e., environmental, social and economic).

More generally, our review showed that the setting up of a more sustainable soil management in viticulture requires a multidisciplinary approach involving all the players throughout the whole production chain. In this context, the availability of plant material more resistant to pathogens and physiologically more suitable for the edaphic environment is essential. Moreover, Cu-induced antibiotic resistance in soil microorganisms highlights the relevance of the *One Health* concept, where the protection of plant, animal and human health is considered to be one. For instance, the frequent contamination of the water sources used for crop irrigation with synthetic compounds for human care (such as drugs like ibuprofen, [152] clearly indicates that the defense plan against pathogens of agricultural crops like vines using agrochemicals should not be considered independently of animal and human health, and vice versa. Actually, thanks to an approach like this, it is possible to effectively contribute to both the continuation of the viticulture (and more generally that of agricultural production) and the preservation of the environment as a whole, in particular the soil ecosystem.

**Author Contributions:** S.C., T.M. and G.O. conceived and designed the review; S.C., Y.P., L.B., G.O., P.L., F.M., G.G., M.S., G.B., R.T., G.V. and T.M. wrote the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by grants from Free University of Bolzano (TN5056, TN2612, TN2081). We acknowledge support by the unibz Open Access publishing fund, Free University of Bolzano.

**Conflicts of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

#### **References**


## *Article* **Recycling Biogas Digestate from Energy Crops: Effects on Soil Properties and Crop Productivity**

**Roberta Pastorelli <sup>1</sup> , Giuseppe Valboa <sup>1</sup> , Alessandra Lagomarsino <sup>1</sup> , Arturo Fabiani <sup>1</sup> , Stefania Simoncini 1,† , Massimo Zaghi <sup>2</sup> and Nadia Vignozzi 1,\***


**Abstract:** Digestate from biogas production can be recycled to the soil as conditioner/fertilizer improving the environmental sustainability of the energy supply chain. In a three-year maize-triticale rotation, we investigated the short-term effects of digestate on soil physical, chemical, and microbiological properties and evaluated its effectiveness in complementing the mineral fertilizers. Digestate soil treatments consisted of combined applications of the whole digestate and its mechanically separated solid fraction. Digestate increased soil total organic C, total N and K contents. Soil bulk density was not affected by treatments, while aggregate stability showed a transient improvement due to digestate treatments. A decrement of the transmission pores proportion and an increment of fissures was observed in digestate treated soils. Soil microbial community was only transiently affected by digestate treatments and no soil contamination from Clostridiaceae-related bacteria were observed. Digestate can significantly impair seed germination when applied at low dilution ratios. Crop yield under digestate treatment was similar to ordinary mineral-based fertilization. Overall, our experiment proved that the agronomic recycling of digestate from biogas production maintained a fair crop yield and soil quality. Digestate was confirmed as a valid resource for sustainable management of soil fertility under energy-crop farming, by combining a good attitude as a fertilizer with the ability to compensate for soil organic C loss.

**Keywords:** anaerobic digestion residues; soil amendment; soil fertilization; soil organic C; soil porosity; soil microbial community

#### **1. Introduction**

Interest in biogas production has grown significantly in the past two decades, following the need to reduce fossil fuel consumption in favour of renewable energy sources. To encourage biogas market penetration, EU policy issued financial incentives [1] which have led to a significant increase in the number of biogas plants. More than 18,000 biogas plants were registered by October 2020 [2] with an overall installed electric capacity (IEC) of 13,520 MW estimated at the end of 2019 [3]. Currently, in Europe, Italy and Germany rank first in terms of the number of active biogas plants, with most Italian plants (1900 units with an IEC of around 1000 MW) located in the Po Valley and other northern regions [4,5].

Biogas production from anaerobic digestion mainly relies on four types of biomass sources: (i) biomass wastes from farms (animal slurries and crop residues) and households (municipal solid waste and food waste); (ii) agro-industrial by-products; (iii) sewage sludges; (iv) biomass from dedicated energy crops [6].

**Citation:** Pastorelli, R.; Valboa, G.; Lagomarsino, A.; Fabiani, A.; Simoncini, S.; Zaghi, M.; Vignozzi, N. Recycling Biogas Digestate from Energy Crops: Effects on Soil Properties and Crop Productivity. *Appl. Sci.* **2021**, *11*, 750. https:// doi.org/10.3390/app11020750

Received: 21 December 2020 Accepted: 11 January 2021 Published: 14 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The energy derived from anaerobic digestion is considered to be almost "carbonneutral" and to bring environmental and social benefits, contributing to a reduction of greenhouse gas emissions (allowed by replacement of fossil fuels) and organic wastes [7], and supporting rural development and new employment opportunities [8]. Against these benefits, biogas production from energy crops generates several issues and conflicts that are under the political attention on a world scale, since the shift of farmland to nonfood systems creates doubts concerning the security of food supply and the environmental impact of energy crops cultivation [9]. One main concern is the environmental sustainability of energy crop cultivation as large amounts of organic matter and plant nutrients are removed with the crop biomass from the field. Depletion of organic matter and plant nutrients from the agricultural system can lead to soil degradation if not balanced by appropriate replenishments. Secondly, since the number of biogas plants in many European countries has increased significantly in recent years, the disposal of residues from anaerobic digestion has become of growing concern [10]. From a sustainability perspective of the biogas supply chain, since a wide range of undecomposed organic compounds and plant nutrients removed from the field (mainly ammonia and phosphate) are retained in the digestate [11–14], the direct land application of digestate is an economical option for residues disposal and soil amendment/fertilization [15–17]. The risk of a potential transfer of organic pollutants, such as herbicides and fungicides, from digestate to rotational crops and feedstuffs is considered very low by the European Food Standards Agency [16].

A third concern is that energy crops require resources (land, water and energy) which inevitably become no longer available for food production [8,9,18,19]. For cereal crop-based productions, the "food vs. fuel" conflict would be overcome if the grain was excluded from the biogas feedstock and used for livestock feed, while in general, the conflict would not exist at all if energy-crop cultivation was carried out on soils unsuitable for food production (marginal land) [20]. In this context, energy crop farming is an effective and profitable strategy to prevent the land from abandonment and degradation while promoting rural investments and new job opportunities [8].

Digestate can be applied to the soil without further processing (whole digestate, WD) [17] or after mechanical separation to obtain a solid fibrous material (solid digestate, SD) which can be directly spread to the field, composted, or dried for intermediate storage and transport [17,21]. Both WD and SD are sources of organic carbon and plant nutrients but since they exhibit quantitative and qualitative differences, they are expected to contribute differently to soil organic matter turnover [22], plant nutrients availability, and soil physical properties [23]. Typically, SD exhibits a great percentage (38–75%) of highly stable organic matter and a low NH4-N to total-N (TN) ratio [23], which make it suitable for use as a soil conditioner rather than as a source of readily available N. The use of digestate as a soil amendment can contribute to soil carbon sequestration, especially in intensively cultivated soils where crop residues are removed [24]. Organic matter addition is beneficial to soil fertility, since it may improve soil structure, increase plant nutrient retention, and water holding capacity and stimulate microbial activity [25]. A higher microbial activity, in turn, may enhance the release of plant nutrients from added residues and soil organic matter itself [26]. Conversely, the low organic matter concentration and the high NH4-N/TN ratio in WD makes it more suitable for use as an N-fertilizer [22,23]. The efficiency of digestate as N fertilizer changes with the features of digestate itself, soil type, crop and time of spreading [4]. Like any other fertilizer, WD should be applied at appropriate rates and times during the crop growing season, to ensure optimum plant nutrient uptake and to avoid phytotoxic effect and pollution of groundwater [16].

Research on digestate suitability for land application is relatively recent and is focused, on the one hand, on agricultural benefits of digestate as soil fertilizer and/or improver, and on the other hand, on the environmental risks associated with digestate use. Overall, many studies have investigated the potential of digestate as N fertilizer and the fate of N in the soil after land application [27,28], as well as the effect of digestate on soil organic matter and chemical properties [28–30], while there is still little knowledge about the

impact of digestate on soil physical [27,28,31] and biological properties (bacterial and fungal communities) [32,33], which are key factors of soil functioning. Knowledge gaps about appropriate rates and soil-digestate interactions still exist, and the research field is very broad and complex, involving different kinds of feedstocks, crops, soils, environments, and agricultural management.

The main goal of our research was to understand the short-term effects of digestate on soil properties through a holistic approach, investigating soil physical, chemical, and microbiological properties and their interactions. Furthermore, we evaluated the effectiveness of digestate in replacing mineral fertilizers and as a resource to compensate for carbon depletion due to biomass removal in a three-year energy crop rotation. The study included both the whole digestate and its mechanically-separated solid fraction.

#### **2. Materials and Methods**

#### *2.1. Digestate*

The whole digestate was obtained from the biogas plant of the Cooperativa Agroenergetica Territoriale (CAT) in Correggio (Reggio Emilia, Italy). The digester was fed with the above-ground biomass from energy crops, including maize, triticale and sorghum silages, combined with by-products from the agricultural industry (i.e., stalks of grapes and sugar beet pulps), and cattle slurry from Parmesan cheese farms [34]. The solid fraction (SD) was retrieved from the whole digestate (WD) through a mechanical solid/liquid separation system following the digestion. SD was rich in organic C (44.4% of air-dry digestate) but relatively poor in N, P and K, whereas WD had a very low organic C content (1.1% of air-dry digestate) and a low C/N ratio (3.1), with a large proportion of NH4-N in the total amount of N (about 60% of air-dry digestate) (Table S1).

For a more in-depth characterization of digestate, we performed molecular-level analyses of microbial communities (see the paragraph on soil sampling and analysis) on WD, SD and two additional fractions, one collected directly from the fermentation silos, the other one from the mechanically-separated liquor (LD).

#### *2.2. Study Site and Experimental Design*

The experimental field was a 35 × 130 m area belonging to the R.G.R. Farm (CAT cooperative partner) located in the lower Po Valley (Correggio, Reggio Emilia, Italy; 44◦49′ N– 10◦45′ E). The land use of the area had been converted from sugar beet cultivation to a 2-year maize-triticale rotation to feed the biogas plant, according to the set-aside scheme introduced by the Common Agriculture Policy. The trial was carried out from January 2011 to October 2013, maize cultivated from spring to summer 2011 and 2013, and triticale from autumn 2011 to summer 2012. The effects of digestate application on soil properties were investigated in the two maize growing seasons, using the whole digestate (WD) as a partial or total replacement of mineral fertilizer, and the digestate solid fraction (SD) as a soil amendment. Nitrogen fertilization was performed during maize post-emergency stage as follows: D0 plots, with mineral fertilizer only (control); D50 plots, based on WD + mineral fertilizer; D100 plots, with WD only. The SD fraction was applied to the WD-fertilized plots (D50 and D100) between one crop cycle and the next.

The treatments were assigned to 4 × 10 m field plots according to a randomized block design with three blocks (replicates). 1 m between plots and 5 m between blocks were kept free to avoid disturbance during soil tillage and to allow machinery operations. During the trial period, the mean annual air temperature was 14.2 ◦C and precipitation 681 mm (Figure S1). The experimental soil was a Hypocalcic Hypovertic Calcisols [35], with a silty-clay texture (Table S2). The main soil physical and chemical characteristics at the start of the trial (time t0) are given in Table S2.

In September 2010, the field was ploughed and harrowed for seedbed preparation. On 2 April 2011, maize (*Zea mais* L., cv. Kalumet) was sown at a density of 7 plants/m<sup>2</sup> and all plots fertilized with urea (125 kg N ha−<sup>1</sup> ). At the plant emergence (20 May 2011), an additional N fertilization was applied as follows: D0 plots, urea (125 kg N ha−<sup>1</sup> ); D50 plots,

urea (62.5 kg N ha−<sup>1</sup> ) plus WD (17,400 L ha−<sup>1</sup> = 62.5 kg N ha−<sup>1</sup> ); D100 plots, WD only (34,700 L ha−<sup>1</sup> = 125 kg N ha−<sup>1</sup> ). WD was spread on the soil surface along the maize rows using mobile equipment (Figure S2) specifically developed by CAT and Cavazzuti Franco (Carpi, Modena, Italy), consisting of a 1 m<sup>3</sup> tanker mounted on a tractor and connected to a boom with 4 trailing hoses, with a 2.80 m working width. The tanker was equipped with a pump-loading apparatus for filling. Maize was irrigated on 26 May, 13 June and 6 July 2011, and harvested at the wax ripeness stage (17 August 2011). On 16 September 2011, the D50 and D100 plots received 40 m<sup>3</sup> ha−<sup>1</sup> SD (equivalent to 10 t ha−<sup>1</sup> ), applied by a solid manure spreader (Vaschieri, Solignano di Castelvetro, Modena, Italy) and incorporated into the soil by ploughing and harrowing. Triticale (x *Triticosecale* Wittm., cultivar Agrano) was sown in November 2011 at a density of 240 kg seeds ha−<sup>1</sup> and fertilized in a single operation in April 2012 by urea only (30 kg ha−<sup>1</sup> ). Due to the high plant density and the lack of suitable equipment for WD application in the standing crop, no WD top-dressing treatment was possible for triticale. The option of a pre-sowing WD treatment was discarded because of the low N use efficiency in the autumn-winter period and the related risk of N leaching [4]. Triticale biomass was harvested on 24 June 2012. In October 2012, the D50 and D100 plots were amended with the SD fraction (40 m<sup>3</sup> ha−<sup>1</sup> ) and the whole field was prepared for maize sowing as previously described (19 April 2013). The trial continued with a maize cycle according to the same practices as for the first experimental year. Due to unfavourable weather conditions (Figure S1), sowing, fertilization and harvesting operations needed to be delayed for about one month, respectively. Maize was harvested at the wax ripeness stage on 3 September 2013.

The whole above-ground biomass yielded at the end of the crop cycles was harvested and used as feedstock for biogas production.

The combination of both SD and WD with the agricultural management (fertilization factor) and sampling data (time factor) were the factors considered for the evaluation of differences in soil physical-chemical and biological characteristics.

#### *2.3. Seed Germination Bioassay*

Extracts of the two digestate fractions (WD and SD) collected from biogas plant at the beginning of experimentation were prepared by adding 25 g digestate to 100 mL of sterile deionized water. The suspensions were shaken for two hours and then centrifuged at 5000 *g* for 30 min. The supernatants from each digestate were used to prepare test solutions with digestate concentrations of 100% (pure), 50%, 25%, 12.5% and 0% (distilled water as control). Petri dishes of 9 cm diameter were prepared, each containing twenty maize seeds placed upon two sheets of Whatman N. 1 filter paper pre-treated with 10 mL of the test solution. The dishes were transferred to a germination chamber under controlled temperature (20 ◦C) in the dark. There were five replicates for each treatment.

The number of seeds germinated in each Petri dish was counted after three days and after one week of incubation, and the germination index (GI) was calculated as a percentage relative to the control [26]. Seedling root elongation was measured after 1 week.

#### *2.4. Crop Yield*

Crop yield just before harvest (in August for maize and in June for triticale) was estimated by collecting biomass at ground level from three randomly selected point in each plot spaced 30 cm from the edges to avoid border effects. In each sampling point, maize was harvested from 1 m in length row sections (including 6–7 plants), while triticale was harvested from 0.5 m<sup>2</sup> areas. After weighing, the biomass was oven-dried at 70 ◦C until constant weight (about 56 h for maize and 48 h for triticale) to determine the dry weight.

#### *2.5. Soil*

#### 2.5.1. Sampling

Soil samples were collected before maize sowing (25 March 2011 = t0; at the beginning of the trial); after maize harvesting (17 November 2011 = t1); before sowing in the second maize cycle (14 April 2013 = t2); at the end of the trial (3 October 2013 = t3).

For soil chemical, biochemical, microbiological and particle size analyses, each plot was sampled by auger to a depth of 20 cm in three selected points, collecting soil cores of 5 cm in diameter. The three cores were then mixed thoroughly providing a single composite sample per plot (3 replicates for each treatment, as a whole). Before chemical and biochemical analyses, the soil was air-dried, ground and sieved through a 2 mm mesh size. The samples for microbiological analyses were stored untreated at −80 ◦C until analysis.

For soil bulk density (BD) and macro-porosity measurements, three undisturbed soil samples were collected from the central part of each plot, at depth increments of 0–10 and 10–20 cm, using a hammer-driven linear sampler. Samples for BD were collected at each sampling time whereas those for macro-porosity analysis were taken only at t<sup>0</sup> and t2.

Soil aggregate stability was determined at t0, t<sup>1</sup> and t<sup>2</sup> on a single composite sample per plot, obtained from three spatially separated sub-samples of soil aggregates collected down to 10 cm depth.

#### 2.5.2. Chemical Analyses

Soil pH was measured potentiometrically in a 1:2.5 soil-water suspension. Soil cation exchange capacity (CEC) and exchangeable base concentrations (Ca, Mg, K and Na) were determined on BaCl<sup>2</sup> triethanolamine (pH 8.2) extracts by flame atomic absorption spectrometry [36]. Soil available Cu, Zn, Fe and Mn were extracted and quantified according to Lindsay and Norvell [37]. Soil total organic carbon (TOC) and total nitrogen (TN) in the bulk soil were measured by dry combustion using a Thermo Flash 2000 CN elemental analyser (Thermo Fisher Scientific, Walthman, MA, USA). The analysis was performed on 20 to 40 mg of soil weighed into Ag-foil capsules and pre-treated with 10% HCl until complete removal of carbonates.

#### 2.5.3. Biochemical Analyses

Microbial biomass carbon and nitrogen (MBC and MBN, respectively) were determined following the fumigation extraction method [38]. Two aliquots from each soil sample were brought to 60% of water holding capacity (WHC), 24 h before the analysis; a first aliquot was immediately extracted with K2SO<sup>4</sup> (0.5 M) and then filtered with Whatman n. 42 filter paper; the second aliquot was fumigated for 24 h at 25 ◦C with CHCl<sup>3</sup> and extracted as the first one. The organic C and N concentration in the extracts was then determined by Thermo Flash 2000 CN elemental analyser (Thermo Fisher Scientific). MBC and MBN were calculated as the difference between the C and N extracted from the fumigated samples and those extracted from non-fumigated samples, respectively.

Soil microbial respiration was determined according to Badalucco et al. [39]. Each sample was incubated at 28 ◦C in a flask sealed with a stopper. The CO<sup>2</sup> developed during incubation was trapped in NaOH solution after 1, 3, 7, 10, 14, 21 and 28 days and then measured by titration with HCl (0.1 M). The cumulative amount of CO<sup>2</sup> produced over 28 days of incubation (MRcum) was regarded as the potentially mineralizable C.

#### 2.5.4. Microbiological Analyses

Soil RNA was extracted using the RNA PowerSoilTM Total Isolation Kit (MoBio, Solano Beach, CA, USA), following the manufacturer's instructions with the minor modification of adding Na-EDTA (0.5 M) to the lyses solution to improve the DNA desorption from clay particles [40]. RNA was eluted in nuclease-free water (Promega, Madison, WI, USA) and then DNA was co-extracted by the RNA PowerSoilTM DNA Elution Accessory Kit (MoBio). The extracted RNA was subsequently subjected to DNase digestion using the

RQ1 RNase-free DNase (Promega) and complementary cDNA was generated by reverse transcription (RT) using the ImProm-IITM Reverse Transcriptase System (Promega).

For Denaturing Gradient Gel Electrophoresis (DGGE) analysis of microbial communities, the extracted DNA and the generated cDNA were amplified using specific primers for bacterial and archaeal 16S rDNA, and for fungal 18S rDNA (Table S3). Amplification and DGGE procedures were carried out following Pastorelli et al. [41] and Lazzaro et al. [42].

Representative bands from archaea and Clostridiaceae-related DGGEs were excised, eluted from gels and screened according to Pastorelli et al. [43]. Selected bands were subjected to direct sequencing by Macrogen Service (Macrogen Ltd., Amsterdam, The Netherlands). The nucleotide sequences collected in this study were deposited in the GenBank database under the accession numbers MF415444-MF415489.

#### 2.5.5. Physical Analyses

Soil texture was determined by the pipette method [44]. Soil bulk density (BD) was measured by the core method according to Blake and Hartge [45].

Soil macro-porosity was determined by the micro-morphometric method [46]. This method allows the characterization and quantification of soil macro-porosity according to pore shape, size distribution, irregularity, orientation and continuity from vertically oriented thin sections of 5.5 × 8.5 cm size, obtained from undisturbed soil samples. A 2.82 × 3.54 cm area of each thin section was captured with a video camera avoiding the edges where disruption could have occurred. The images collected were then analysed by Image-Pro Plus software (Media Cybernetics, Silver Spring, MD, USA), set up specifically to measure pores larger than 50 µm. The total porosity and pore distribution were calculated from the measurement of pore shape and size [46]. From a functional point of view, the elongated pores of 50–500 µm were described as transmission pores and the pores with >500 µm size as fissures [47]. The thin sections were also examined by a Zeiss "R POL" microscope at a 25× magnification to characterize soil structure.

Soil aggregate stability was determined by the wet sieving method and the calculation of the mean weight diameter of water-stable aggregates (MWD) [48]. Soil aggregates from each composite sample were air-dried, weighed and separated into different size fractions (10–20, 4.75–10, 2–4.75, 1–2, <1 mm) using a vibrating sieve shaker (Retsch, Germany). The most representative aggregate size fraction was used to perform wet sieving. Twenty grams of aggregates from the most abundant size class (4.75–10 mm) were directly soaked for 5 min on the top of nests of 4.75, 2, 0.25- and 0.05-mm diameter sieves immersed in water (fast wetting). The nest of sieves with its content was then vertically shaken in water by an electronically controlled machine with a stroke of 40 mm per 10 min, at a rate of 30 complete oscillations per minute. For each sample, 3 repetitions were performed.

#### *2.6. Statistical Analyses*

The results of soil physical, chemical and microbiological (richness and α-diversity indices) analyses were processed by analysis of variance (ANOVA) followed by Fisher's least significant difference (LSD) test at the significance level *p* ≤ 0.05, using the Statistica software (Palo Alto, CA, USA). Pearson correlation analysis was performed among physicochemical properties of soil by Statistica software.

Band migration distance and intensity for each DGGE profile were obtained using the Gel Compare II software v 4.6 (Applied Maths, Saint-Martens-Latem, Belgium). The number of bands (species richness) and their relative abundance (Shannon index, H' and Simpson index, D) were used as proxies of richness and α-diversity of soil microbial communities, as described by Pastorelli et al. [43]. The banding patterns of bacterial and fungal DGGE profiles were converted into presence/absence band matching tables and imported into PAST3 software [49]. Non-metric multidimensional scaling (MDS) based on the Dice coefficient was performed to represent the distance between the DGGE profiles in the two-dimensional space. Analysis of similarity (ANOSIM) based on Dice similarity coefficient and 9999 permutational tests were run to assess the statistical significance in microbial community structure due to fertilizer/amendment treatments.

Nucleotide sequence chromatograms were edited using Chromas Lite software v2.1.1 (Technelysium Pty Ltd., Tewantin, Old, AU) to verify the absence of ambiguous peaks and to convert them to FASTA format. The DECIPHER's Find Chimeras web tool [50] was used to uncover chimaeras hidden in nucleotide sequences. The Web-based BLAST tools was used to identify closely related nucleotide sequences within those stored in the GenBank database.in order Microbial taxonomic identification was achieved by means of different sequence similarity thresholds as described by Webster et al. [51].

#### **3. Results**

#### *3.1. Germination Index Bioassay*

GI was lowest when maize seeds were treated with undiluted SD and WD soluble extracts (57% and 34.9%, respectively; Table 1) but increased with increasing dilution ratio. According to McLachlan et al. [52], GI values above 70%, as those observed under 50%, 25% and 12.5% digestate concentrations, indicate the absence of toxicity.

**Table 1.** Relative seed germination index (GI) of maize under different digestate concentrations, expressed as the percentage of germinated seeds relative to the control GI (distilled water).


SD = solid digestate; WD = whole digestate.

After one-week incubation there were significant differences in rootlet lengthening between SD and WD treated maize seeds, as well as between the different digestate used concentrations. The rootlet length was lowest in the undiluted extracts and greatest under 12.5% (SD) and 25% (WD) concentrations (Table 2). With a 12.5% SD concentration, the root elongated more than in the control, although not significantly.

**Table 2.** Root length (mm) of maize seedlings after 1 week of incubation under different digestate concentrations (means followed by standard error in brackets). Different Latin letters within a column indicate statistically significant differences between digestate concentrations at *p* ≤ 0.05 (Fisher LSD test); different Greek letters within a row indicate statistically significant differences between digestate types at *p* ≤ 0.05 (Fisher LSD test).


SD = solid digestate; WD = whole digestate.

#### *3.2. Crop Biomass Yield*

Neither maize nor triticale biomass showed significant differences between treatments (Table 3). In 2013, due to abundant rainfall (Figure S1), the growth of maize suffered a marked delay compared to 2011, along with a reduction in the biomass yield irrespective of treatment (8.8–12.8 t ha−<sup>1</sup> against 19.6–22.1 t ha−<sup>1</sup> , respectively).


**Table 3.** Maize and triticale above-ground biomass (t dry weight ha−<sup>1</sup> ) under different experimental treatments (means followed by standard error in brackets).

D0 = 100% N as urea; D50 = 50% N as urea + 50% N as WD; D100 = 100% N as WD.

#### *3.3. Soil Chemical Properties*

As reported in Table 4, the average TOC content in the plots under digestate treatment generally showed slight increases compared to that of the control plots, with no difference between the application rates.

**Table 4.** Soil total organic carbon (TOC), total nitrogen (TN) and C/N ratio under the different experimental treatments at different sampling times (means followed by standard error in brackets). Different letters indicate statistically significant differences at *p* ≤ 0.05 (Fisher LSD test).


D0 = 100% N as urea; D50 = 50% N as urea + 50% N as WD; D100 = 100% N as WD.

TN followed a different trend but, overall, it was well correlated with TOC (0.723 \*\*\*), confirming a positive digestate effect in the third experimental year (t<sup>2</sup> and t<sup>3</sup> sampling). The TOC to TN ratio did not change with treatment, except for t<sup>1</sup> sampling which showed lower C/N values under digestate application (Table 4).

Soil CEC values, exchangeable Ca and exchangeable Mg concentrations did not differ among treatments for the entire duration of the trial (Table 5). In contrast, at the end of the first year (t1) the exchangeable K concentration was increased by D50 and D100 regardless of the application rate. At t2, K showed a significant decrease in all plots as compared to t<sup>0</sup> and t<sup>1</sup> contents. The available Cu, Zn, Fe and Mn contents were not affected by treatments (Table 5).

**Table 5.** Soil cation exchange capacity (CEC), exchangeable bases (K, Na, Mg, Ca) and available metal content (Cu, Zn, Fe, Mn) in the soil (depth = 0–20 cm) under different fertilization treatments (D0, D50, D100) and at different sampling times (t<sup>0</sup>, t1, <sup>t</sup>2) (means followed by standard error in brackets). Different letters indicate significant differences between soil samples at*p*≤0.05 (Fisher LSD test).


D0 = 100% N as urea; D50 = 50% N as urea + 50% N as WD; D100 = 100% N as WD.

#### *3.4. Soil Physical Properties*

Soil BD did not change significantly with treatments and was stable over time (Figure 1).

**Figure 1.** Soil bulk density (BD; g cm−<sup>3</sup> ) at 0–10 and 10–20 cm depth, under the different treatments and at the different sampling times.

Macro-porosity ranged within moderate levels in the surface layer (10–25%) while it averaged less than 5% in the deeper layer, indicating a very compact soil, as defined according to the micro-morphometric method [53] (Figure 2). Differences between treatments were significant in the surface layer only. The t<sup>0</sup> sampling showed a certain degree of field variability for soil macro-porosity, with D0 plots showing a higher macro-porosity than D50 and D100 plots (related to a larger number of fissures) and D50 plots featuring a higher proportion of irregular pores compared to D0 and D100 plots. In the t0–<sup>2</sup> time frame, soil total macro-porosity increased under D100 with an increase in the percentage of >500 µm elongated pores (fissures) and a reduction in that of 50–500 µm elongated pores. Over the same period, macro-porosity remained quantitatively unchanged under D50, showing a decrease in the 50–500 µm elongated pores.

**Figure 2.** Soil macroporosity (pores size >50 µm) expressed as a percentage of area occupied by pores of different shape (regular, irregular and elongated pores) at 0–10 cm (**A**) and 10–20 cm (**B**) depth and at two different sampling times (t<sup>0</sup> and t<sup>2</sup> ). Different letters above bars indicate statistically significant differences between the % of fissures (elongated pores, size > 500 µm) in relation to the total macro-porosity; different letters inside the bars indicate significant differences within each shape or size class of pores at *p* ≤ 0.05 (Fisher LSD test).

Soil aggregate stability was very low at the beginning of the trial (MWD < 2.5 mm, against a theoretical MWD maximum of 7.375 mm for the 4.75–10 mm size class aggregate)

but increased over time regardless of treatment (Figure 3). The effect of digestate treatment was significant only at t1, soon after WD distribution. After two years (t2), the differences in soil aggregate stability between treatments were not significant.

**Figure 3.** Soil aggregate stability as expressed by the mean weight diameter index (MWD; mm), under the different treatments and at the different sampling times. Different letters indicate statistically significant differences between treatments and sampling times at *p* ≤ 0.05 (Fisher LSD test).

#### *3.5. Soil Microbial Biomass and Respiration*

Soil MBC was relatively stable over time with a small but significant increase (*p* ≤ 0.05) only in D50 plots. The MBN decreased from t<sup>1</sup> to t<sup>2</sup> regardless of treatments (Table 6).

**Table 6.** Soil microbial biomass C (MCB), microbial N (MBN) and cumulative microbial respiration (MRcum) under different fertilization treatments (D0, D50, D100) and at different sampling times (t<sup>0</sup> , t1 , t<sup>2</sup> ) (means followed by standard error in brackets). Different letters indicate significant differences between soil samples at *p* ≤ 0.05 (Fisher LSD test).


D0 = 100% N as urea; D50 = 50% N as urea + 50% N as WD; D100 = 100% N as WD.

The C mineralization potential (after 28 days of incubation) did not change significantly either in relation to treatment or time, except for D50 plots where it was higher than in the control plots at t<sup>1</sup> (Table 6).

#### *3.6. DGGE Analysis of Total Bacterial and Fungal Communities*

The abundance (richness) and α-diversity (Shannon–Weiner and Simpson indices) calculated from DGGE profiles showed that the soil bacterial community was overall richer and more diverse than the fungal community (Table S4). When considering all groups independently (12 groups: 4 sampling time combined with three digestate treatments), there were significant differences between soil samples for both bacterial and fungal communities (Table S4). Multifactorial ANOVA (Table S4) showed that the species richness and α-diversity indices of the bacterial community were significantly influenced by the interaction between sampling time and digestate treatment. Differently, only the sampling time had a significant effect on the species richness and α-diversity indices of the fungal community (Table S4).

At t0, MDS ordination showed a low inter-specific variation between the bacterial communities from the differently treated plots (Figure 4A). At t1, the D50 and D100

bacterial communities were clearly separated from the D0 ones, which grouped with t<sup>0</sup> communities. Bacterial communities at t<sup>2</sup> and t<sup>3</sup> grouped together regardless of treatment and were well separated from the t<sup>0</sup> and t<sup>1</sup> ones (Figure 4A). Conversely, at t<sup>0</sup> the fungal community showed a higher inter-specific variation than bacterial community. In the following sampling, the fungal community showed a progressive change of its structure, which seems to be independent of the treatments (Figure 4B).

**Figure 4.** MDS ordination plots of bacterial 16S rDNA (**A**; stress = 0.218) and fungal 18S rDNA (**B**; stress = 0.282). Symbols: circle = D0 treatment; triangle = D50 treatment; square = D100 treatment. Colours: white = time t<sup>0</sup> ; red = time t<sup>1</sup> ; blue = time t<sup>2</sup> ; black = time t<sup>3</sup> .

Due to the poor reliability of MDS ordination results, especially for 18S-DGGE (stress = 0.282), DGGE profiles were further analysed by multivariate analysis. When testing all groups independently (sampling time × fertilizer treatment), the one-way ANOSIM global test revealed significant differences in both bacterial and fungal DGGE profiles (Table 7), although R values were not sufficiently reliable. According to the outcomes of the two-way crossed ANOSIM test, the differences in both bacterial and fungal community composition were greater in relation to the sampling time (R = 0.822 and 0.808 for bacteria and fungi, respectively) than in relation to the digestate treatment (R = 0.448 and 0.275 for bacteria and fungi, respectively) (Table 7).


**Table 7.** Summary of ANOSIM analysis based on 16S- and 18S-rDNA Dice similarity matrices. In the one-way ANOSIM groups were analysed independently (three digestate treatments vs. four sampling time), whereas the two factors (sampling time and digestate treatments) were analysed by a two-way analysis.

#### *3.7. DGGE Analysis of Total Active Bacterial Community*

The active bacterial community was analysed by matching the t<sup>0</sup> DGGE profiles with those obtained at t<sup>1</sup> (different seasons within the same year: March vs. November) and t<sup>2</sup> (different years under the same field conditions: before maize sowing, 2011 vs. 2013).

The abundance (richness), α-diversity (Shannon–Weiner and Simpson indices), and composition of the active bacterial community were more influenced by sampling time than by digestate treatment or sampling time × digestate treatment interaction (Data not shown). The separation between active bacterial communities was stronger when they were compared according to the different season (t<sup>0</sup> vs. t1) than to the different year (t<sup>0</sup> vs. t2) (Table 8).

**Table 8.** Summary of ANOSIM analysis based on the 16S-cDNA Dice similarity matrices. In the one-way ANOSIM groups were analysed independently (three fertilizer treatments vs. four sampling time), whereas the two factors (sampling time and fertilizer treatments) were analysed by a two-way analysis.


#### *3.8. DGGE Analysis of Archaea and Clostridiaceae-Related Communities*

The DGGE profiles from the different digestate fractions were very similar to each other and quite different from those of the soil (Figure S3). Digestate-based treatments had no substantial effect on soil archaeal (Figure S3a) and Clostridiaceae-related bacterial (Figure S3b–e) communities. Some additional dominant bands were found at t<sup>1</sup> in D50 and D100 DGGE profiles obtained with the primer sets specific for Clostridiaceae-cluster I and -cluster IV. In particular, a group of γ-Proteobacteria-related bands appeared in the Clostridiaceae-cluster I DGGE profiles (Figure S3b), while a group of bands phylogenetically related to *Caproiciproducens galactitolivorans* (similarity ranged from 93% to 94%) appeared in the Clostridiaceae-cluster IV DGGE profiles (Figure S3d). These bands were almost undetectable in the t<sup>2</sup> DGGE profiles.

Overall, the digestate DGGE profiles were characterized by one or more all-time dominant bands related to Clostridiacea (Figure S3b–e; Table S5). However, none of the primer sets was specific enough to detect only *Clostridium*-related species, since several DGGE bands revealed to be related to β-, δ- and γ-Proteobacteria divisions, Acidobacteria group or Actinobacteria phylum (Table S5).

#### **4. Discussion**

#### *4.1. Effects of Digestate on Soil Chemical, Physical and Microbiological Properties*

In this trial, digestate treatments provided consistent results in the two years of maize cultivation. The soil TOC tended to be slightly higher in plots treated with digestate than in plots under mineral fertilization, in agreement with results obtained using digestate or other different bioenergy by-products as a soil amendment or fertilizer [28,29]. It is possible that, to some extent, soil organic C enrichment was limited by tillage practices, due the dilution of the organic matter across the ploughed layer and the exposure of physically protected organic compounds to enhanced mineralization [54].

Functional properties of organic residues as amendments are related to the organic matter stability, i.e., the ratio of recalcitrant to labile organic fractions [55] and how these interact with soil features, climate and crop management. There is consistent laboratory evidence of lower carbon mineralization of digestate compared to undigested feedstock, due to an increase of the recalcitrance of organic matter during digestion [28]. However, results from previous short-term experiments on the effects of digestate on soil carbon and nitrogen and crop yield are contrasting, probably due to the various chemical characteristics of digestate and different type of soil used in the experimentations [13,33,56,57].

Overall, the role of soil organic matter in soil fertility and plant nutrition may be summed up in its ability to supply and store plant nutrients [58]. This role is expressed through the release of organically-bound plant nutrients by microbial mineralization and the contribution of organo-mineral complexes to the retention of plant nutrients as available cations [58]. As indicated by the close similarity between soil TN and TOC distribution patterns, soil organic matter contributed to the overall N pool. In addition, the determination coefficient of the relationship between soil TN and TOC (R<sup>2</sup> = 0.723 \*\*\*) suggests that additional factors may account for TN variations, namely the dynamics of mineral N supplied by fertilizers (WD, urea) and soil organic matter mineralization.

According to crop yield performances, digestate treatments were at least as effective as mineral fertilization in supporting crop requirements. There was no evidence of a significant contribution of the organic fraction of the digestate to the cation exchange capacity of the soil, which is explained by the modest TOC variations found after digestate treatments and the fine-textured composition of the soil mineral fraction [59]. This agrees with previous studies showing that the effects of organic amendment on soil CEC were generally more pronounced in coarse-textured soils than in clayey soils [59].

The whole digestate (WD) also proved to be a valuable source of K, by increasing the available K pool of the soil by 22% during the first year of trial. To further support the high potential of digestate as a substitute for mineral K-fertilizers, numerous experimental evidences demonstrate very high rates of K recovery during anaerobic digestion (above 94%) from a wide range of feedstock materials [60].

The unexpected decrease in soil K and Na content across the experimental field before sowing in the second maize cycle (compared to their average content in the previous sampling times) can be ascribed to a leaching effect (Figure S1), which conversely left Ca and Mg concentrations unchanged due to their lower water-solubility and the buffering effect of soil carbonates [61].

Soil BD was not affected by digestate treatments, which disagrees with results by other authors who found a reduction of BD under organic amendments in both compacted and uncompacted soils [62]. BD and organic matter are linked by a close relationship involving physical and chemical interactions between organic substances and soil mineral particles [63]. Usually, due to a lower density of the organic matter compared to that of the mineral fraction, the average BD of a mixture mineral fraction/organic matter decreases as the organic matter content increases. In the present experiment, several factors may have interfered with these relationships, i.e., an experimental period too short compared to the time required for soil structure formation and a contrasting effect of soil tillage on aggregates formation and stabilization. This was reflected in the pattern of soil pore size distribution, with a decrease in the amount of transmission pores, which are of primary importance for optimal soil–water–plant relationships, and an increase in the proportion of fissures mainly involved in water and air flows but related to poor structure and physical degradation when they are (as in D100) over 70% of total porosity [64].

The stability of soil aggregates is a key indicator of soil physical quality, affecting the ability of the soil to retain its structure and the related physical and hydraulic functions against degradative forces [65]. Soil aggregate stability relies on a complex range of factors involving soil texture and mineralogy, the chemistry of soil cations and soil organic matter content and quality [64]. At the beginning of the trial, aggregate stability (expressed as MWD) was quite low possibly due to the high silt proportion and the low organic C content [66,67]. However, aggregate stability was improved by digestate treatment during the first experimental year, consistently with findings of other authors [27,31]. In addition, it correlated positively with TOC (Figure 5), in line with the role of soil organic C as a driver of soil aggregate formation [66–68].

With respect to soil biological parameters, biochemical analysis revealed just a slight (statistically not significant) increase over time in the soil MBC under digestate treatment. This increase was consistent with the trend of TOC, suggesting that part of the organic C supplied by digestate could have been converted into MBC [69]. The small extent of MBC increase was expected from a short term digestate treatment, due to the relatively high recalcitrance of the organic matter in SD and the low organic C concentration in WD [26,70]. This evidence confirms SD application as a valuable tool to improve soil C sequestration [71] and to compensate for C depletion associated to crop biomass removal.

**Figure 5.** Correlation between soil aggregate stability based on the mean weight diameter (MWD) and total organic carbon (TOC) (\*\* *p* ≤ 0.01).

Microorganisms are crucial for soil fertility [25]. They drive the turnover of organic substrates and their abundance and diversity can be affected by soil management as well as by the quality of soil amendments/fertilizers. A reduction in microorganism number and diversity may impair their ability to perform specific functions as well as to withstand soil perturbations from a long-term perspective [72]. In this trial, "time" was the main factor affecting the α-diversity of soil microbial communities. Differences in the microbial community structure in response to digestate addition were showed when t<sup>0</sup> and t<sup>1</sup> were compared. Conversely, the microbial community structure remained quite similar when t<sup>2</sup> and t<sup>3</sup> were compared. In addition, active bacterial communities resulted more affected by season than by digestate treatments, contrasting many reports which indicated an enhanced soil microbial activity after field applications of digestates [32,33]. Calbrix et al. [73], in a study dealing with the impact of organic amendments on soil bacterial communities over a 12-month period, observed that changes in soil bacterial community structure were only temporary and that seasonal variations had the greatest effect on microbial community composition. Accordingly, in our study, digestate showed to have only a transient effect on the microbial community structure. Successively, the soil microbial communities developed new stability and equilibrium over time in both digestate treatments, thus strengthening the hypothesis of a resilience of microbial communities to anthropogenic changes [74–76]. Likewise, the Archaea and Clostridiaceae-related bacteria revealed remarkably stable soil resident communities, with negligible and temporary changes after the introduction of allochthones species (Figure S3).

#### *4.2. Effects of Digestate on Seed Germination and Crop Yield*

The GI bioassays revealed that highly concentrated SD and WD extracts impaired seed germination, whereas <50% digestate concentrations had no phytotoxicity. This suggests that the use of digestate should follow appropriate rates and timings of application to avoid the direct contact with seeds, as also described by Alburquerque et al. [26]. According to our experimental plan, we can exclude any phytotoxic interference of digestate with maize seed germination under field condition, as the SD fraction was applied several months before sowing and WD in the post-emergence stage.

Interestingly, the 12.5% SD concentration increased the germination index and the relative root elongation as compared to the control, which can be explained by assuming a stimulating effect of plant nutrients, growth enhancers or even phytormone-like compounds contained in SD as suggested by other authors [26,77].

In the first two years of the trial, both maize and triticale biomass yields were consistent with the average yields in the area [78], which is promising in the perspective of agricultural use of digestate, alone or combined with mineral fertilizers. The implementation of conservation tillage management [79] may further improve the efficiency of digestate as an amendment and/or fertilizer.

With regard to the pronounced decrease in maize yield across the whole experimental field in the third experimental year (second maize growth season), it was most likely due to a combination of adverse climate and soil physical conditions arising from the abundant rainfall between January and April (Figure S1), which caused a shift of the entire growing season. At the same time, the fine texture of the soil combined with the low organic matter content might have led to a stronger impact of the heavy machinery on soil structure and hydrological behaviour, resulting in insufficient drainage, extended water stagnation and overall poor soil physical conditions for seed germination and shoot development [80].

#### **5. Conclusions**

With a focus on the environmental sustainability of the bioenergy supply chain, the application of digestate to the soil can meet the need to safely dispose and recycle the residues coming from anaerobic digestion and, at the same time, to compensate for soil C and plant nutrient depletion due to crop biomass removal.

From our results, digestate application in a three-year maize-triticale rotation cycle proved to be as effective as 100% mineral fertilization in maintaining crop productivity level. Moreover, the increase in soil TOC following digestate treatments confirmed digestate effectiveness to compensate for carbon depletion.

Further research is needed to increase the knowledge about the optimum dose of digestate to be applied in relation to soil/crop specificities and the best application method to minimize potential negative effects of digestate to the soil and environment quality. The pattern and extent with which the effects of digestate treatments were expressed and their temporal fluctuations underline complex dynamics of chemical, physical and biological processes affecting the material brought to the soil. This suggests that a more or less long period of time is needed during which the achievement of a new stable equilibrium in the soil functions is regulated by the interaction between the amount and quality of biomass supplied, the impact of mechanical operations associated with crop management and climate trend.

Further expected benefits from digestate as amendment, such as improvement of soil bulk density and porosity, were not observed, possibly due to a counteracting interference of soil tillage operations. The effectiveness of a soil amendment and the sustainability of the use of digestate can be strongly conditioned by the crop management system as a whole, and in particular, by those cultivation practices that have a direct impact on the soil and the dynamics of the organic matter and nutrients supplied with the treatment. For this reason, to fully exploit digestate potential, its use should be integrated within an overall more conservative soil management system, involving reduced soil mechanical disturbance. This would be essential to prevent soil physical degradation and excessive organic matter mineralization, thus allowing the organic compounds of digestate to perform their chemical, physical and biological functions and minimize the risk of N loss by leaching and/or gas emissions.

Additional considerations regard the cultivation of energy crops in marginal lands or set-aside areas; this could be a solution to the "food vs. fuel conflict" and, at the same time, would promote rural investments and new job opportunities.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/2076-341 7/11/2/750/s1. References [81–94] are cited in the supplementary materials.

**Author Contributions:** Conceptualization, R.P., G.V., A.L., A.F. and N.V.; formal analysis, R.P., G.V., A.L., S.S. and N.V.; investigation, R.P., G.V., A.L. and N.V.; resources, M.Z.; data curation, R.P., G.V., A.L. and N.V.; writing—original draft preparation, R.P., G.V., A.L. and N.V.; writing—review and editing, R.P., G.V., A.L. and N.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by the Italian Ministry of Agriculture (BIOMASSVAL project).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data is contained within the article or supplementary material.

**Acknowledgments:** We acknowledge the technical support of Azienda Agraria R.G.R. that hosts the field experiment and the CAT cooperative that supplies the digestate. In particular, we greatly acknowledge Gabriele Santi, president of CAT, that coordinated the agronomic management of the field trial. The technical work during sampling of Raimondo Piccolo and Ilaria Secciani is highly appreciated.

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

#### **References**


## *Article* **Partial Least Squares Improved Multivariate Adaptive Regression Splines for Visible and Near-Infrared-Based Soil Organic Matter Estimation Considering Spatial Heterogeneity**

**Xiaomi Wang <sup>1</sup> , Can Yang <sup>2</sup> and Mengjie Zhou 1,\***


**Abstract:** Under the influence of complex environmental conditions, the spatial heterogeneity of soil organic matter (SOM) is inevitable, and the relationship between SOM and visible and near-infrared (VNIR) spectra has the potential to be nonlinear. However, conventional VNIR-based methods for soil organic matter estimation cannot simultaneously consider the potential nonlinear relationship between the explanatory variables and predictors and the spatial heterogeneity of the relationship. Thus, the regional application of existing VNIR spectra-based SOM estimation methods is limited. This study combines the proposed partial least squares–based multivariate adaptive regression spline (PLS–MARS) method and a regional multi-variable associate rule mining and Rank–Kennard-Stone method (MVARC-R-KS) to construct a nonlinear prediction model to realize local optimality considering spatial heterogeneity. First, the MVARC-R-KS method is utilized to select representative samples and alleviate the sample global underrepresentation caused by spatial heterogeneity. Second, the PLS–MARS method is proposed to construct a nonlinear VNIR spectra-based estimation model with local optimization based on selected representative samples. PLS–MARS combined with the MVARC-R-KS method is illustrated and validated through a case study of Jianghan Plain in Hubei Province, China. Results showed that the proposed method far outweighs some available methods in terms of accuracy and robustness, suggesting the reliability of the proposed prediction model.

**Keywords:** soil organic matter; near-infrared spectroscopy; spatial heterogeneity; multivariate adaptive regression splines; partial least squares regression

#### **1. Introduction**

Soil organic matter (SOM) content is significantly relevant to soil fertility, biological productivity, and agricultural sustainable development [1,2]. Accurate prediction of SOM content is of great significance for land management [1,3]. Extensive studies have fully affirmed the capability of SOM prediction methods based on visible and near-infrared (VNIR) spectra [4,5]. However, the function relation between SOM content and VNIR spectra should be established because of the strong soil spatial heterogeneity under the influence of complex environmental conditions. In addition, the formation, variation, and decomposition of SOM are generally influenced by various factors, and the VNIR spectra are comprehensively related to various soil properties; hence, the relation of SOM content to VNIR spectra is complex with high nonlinearity [6,7]. Therefore, further exploration of VNIR spectra-based SOM prediction models is needed to improve simulation and prediction accuracies.

An extensive literature review demonstrates that previous methods for predicting soil properties can be categorized into linear and nonlinear prediction methods. For example,

**Citation:** Wang, X.; Yang, C.; Zhou, M. Partial Least Squares Improved Multivariate Adaptive Regression Splines for Visible and Near-Infrared-Based Soil Organic Matter Estimation Considering Spatial Heterogeneity. *Appl. Sci.* **2021**, *11*, 566. https://doi.org/10.3390/ app11020566

Received: 10 October 2020 Accepted: 4 January 2021 Published: 8 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

multiple linear regression (MLR), partial least squares (PLS) regression, and geographic weighted regression (GWR) are typical linear prediction methods. These methods have been extensively used in many previous studies because they are easily accessible in most software packages [4]. A common characteristic of these prediction methods is that they assume, either explicitly or implicitly, a linear relationship among the analyzed data sets. A nonlinear relationship, however, is prevalent in practice between SOM content and VNIR spectra. Thus, a series of nonlinear prediction methods has been proposed to estimate nonlinear relationships, including, but not limited to, local weighted regression (LWR) [8,9], artificial neural network [10,11], and support vector machine (SVM) methods [12,13]. Although these methods can efficiently determine nonlinear relationships in certain situations, they lead to models based on global optimization and disregard the spatial heterogeneity of the nonlinear relationship. This underlying principle will cause the accuracy of existing methods to easily change with different samples and have difficulty for application to other regions. Therefore, nonlinear estimation methods considering spatial heterogeneity are needed for SOM content prediction.

The multivariate adaptive regression spline (MARS) method is a nonlinear prediction method that considers the effects of local difference. The MARS method is the expansion of splines, thereby providing immense flexibility by automatically determining the splines (i.e., the number of basis functions (BFs) and locations of knots) without human interference. This method has been successfully applied to various fields, such as geotechnical engineering and soil liquefaction assessment [14–16]. The MARS method has proven to be advantageous due to its adaptation to nonlinear relationships and to its adaptive model construction process, while considering the natural local difference in the training data sets. These characteristics of the MARS method indicate its considerable potential to represent the relationships between VNIR spectra and SOM content. Although tracked records of the application of the MARS model to soil property prediction are available in literature [11,17], these records are not well received. The possible reason may be that most of the applications implicitly or explicitly consider the explanatory variables as mutually independent, which is often not the case in reality. In fact, the VNIR spectra of soil samples are the comprehensive representation of soil properties, and thus multicollinearity inevitably exists. As such, multicollinearity should be removed to obtain relatively independent variables when constructing a MARS model. The multicollinearity-removing strategy from PLS regression is verified to be beneficial by reducing the dimensionality, considering the high dimensionality and correlated representations [4,18]. Therefore, this research proposes the utilization of partial least squares regression to remove multicollinearity and obtain extensive information on the explanatory variables; the resultant principal components (PCs) from the PLS regression are used as a replacement for the original explanatory variables to construct a MARS model, thereby overcoming the weakness of the available MARS model for prediction of soil properties. For simplicity, the entire process is referred to as the partial least squares–based multivariate adaptive regression spline (PLS–MARS) method herein. The proposed PLS–MARS method can effectively represent the nonlinear relationship between the data with local difference and numerous variables, as will be discussed in Section 3.

The accuracy of the prediction model relies significantly on calibrated samples. Underrepresentation of calibrated samples will lead to a biased prediction model due to soil spatial heterogeneity. However, preceding studies (e.g., [8,9,11,12]) focused exclusively on developing new prediction models for estimating soil properties, disregarding the influence of the calibration sets on the reliability of the prediction models. In particular, the calibration sets in most of these studies were selected randomly or from several available methods, such as the concentration gradient method (C method) [19,20]. The prediction model may be biased if inappropriate calibration sets are adopted [3]. Hence, reasonable calibration sets should be selected before they are applied to construct prediction models. Accordingly, our newly developed multi-variable analysis method (i.e., regional multivariable associate rule mining and Rank–Kennard-Stone (MVARC-R-KS) method) [3] was

adopted to select calibration sets in the current study (see Subsection III-A for elaboration). The MVARC-R-KS method was used due to the following two reasons. (1) This method is effective in enhancing the accuracy of linear prediction methods, such as PLS regression (Wang, Chen, Guo and Liu [3]), thereby providing the confidence to extend it for nonlinear prediction methods (e.g., the proposed PLS–MARS method in this study). (2) The MVARC-R-KS method can consider multiple influential factors (e.g., chemical component, spectrographic information, and environmental factors) and select a representative calibration set, which makes the VNIR-spectrum–based SOM prediction model calibrated from such a calibration set as to be an extensively applicable one [3]. Three commonly used calibration set selection methods, namely, component concentration representative methods (e.g., the C method [19,20]), spectrum representative methods (e.g., the Kennard–Stone (KS) method [21]), and component concentration and spectrum representative methods (e.g., the Rank–Kennard-Stone (Rank–KS) method [22]), are also utilized and compared with the MVARC-R-KS method to investigate their influences on the prediction models and further verify the effectiveness of the proposed method. To the best of our knowledge, such a systematic comparison of the influences of different calibration set selection methods on the performance of SOM prediction models, particularly nonlinear prediction models, such as the proposed PLS–MARS model herein, appears to be original. This comparison is expected to guide the calibration set selection for a specific prediction model in the future. – – – – – –

–

The remainder of the paper is structured as follows. The study area and materials for this study are introduced in Section 2. The calibration set selection methods and the proposed PLS–MARS prediction method are described in Section 3. Experiments on the simulated data set and real application are discussed in Section 4, illustrating and validating the proposed method. The principal contributions and observations of this study are discussed in Section 5. Section 6 concludes this research. –

#### **2. Study Area and Materials**

#### *2.1. Study Area*

Sampling was done in Jianghan Plain in Hubei Province, China. Jianghan Plain is an important agricultural region because of the humid climate and fertile soil. However, in recent years, the ecological system of the plain has become unstable, and the SOM content has decreased due to long-term agricultural activities. Therefore, estimating the SOM content with high accuracy in this plain is necessary and beneficial to land resource management. To this end, 260 topsoil samples (i.e., 0 cm to 30 cm; see Figure 1) were obtained from this area in June 2014. The samples were taken at least 100 m apart from one another. All samples were partitioned into two parts for chemical study and spectral measurement.

**Figure 1.** Distribution of the soil samples.

#### *2.2. SOM Content and Spectral Measurement*

SOM content was measured by following the analysis and calculation methods in literature [3,23]. An ASD FieldSpec3 portable spectral radiometer was used to measure spectral reflectance, which ranged from 350 nm to 2500 nm with a sampling interval from 1.4 nm to 2 nm. The measurement procedures were carried out in the dark to decrease the interference of external light. The light source was from a standardized white halogen light, which had a 45 ◦ incident angle and was at a distance of 45 cm.

#### **3. Methods**

VNIR spectra-based SOM content prediction generally consists of the following three parts. Part 1 is the spectral preprocessing. Spectral reflectance is unavoidably impacted by random noises, baseline drift, and scattering effects because of the influence of error from the experimental instruments and ambient noises; thus, the stability of the constructed VNIR-spectrum–based prediction model is affected [3]. Hence, spectral preprocessing is a prerequisite for the construction of the prediction model, which will be detailed in Section 3.1. In part 2, the representative calibration set is selected using methods such as the MVARC-R-KS method [3], which will be detailed in Section 3.2. Part 3 aims to build an SOM prediction model based on the VNIR spectra, such as by using typical prediction methods and the proposed PLS–MARS method, which will be elaborated in Section 3.3.1. Lastly, in Section 3.3.2, the accuracy of the constructed prediction model is evaluated by several assessment indices, as will be described in Section 3.3.2. – –

#### *3.1. Spectral Preprocessing*

According to the characteristics of the spectra, spectral preprocessing was conducted as follows. The spectra ranging from 400 mm to 2350 mm was retained to reduce the interference of noises (Liu et al., 2014b; Liu et al., 2016) (Figure 2a). In addition, continuum removed spectral curves exhibited several diminutive absorption valleys (Figure 2b), which could interfere with the prediction based on the VNIR spectra. Hence, Savitzky–Golay (SG) smoothing was used to denoise this interference. Furthermore, multiplicative scatter correction (MSC) and mean center (MC) operations were utilized to reduce the influence of the unavoidable scattering. In short, the spectra were preprocessed successively through SG smoothing, MSC operation, and MC operation. –

**Figure 2.** (**a**) Pretreatment spectral curves of the soil samples. (**b**) Continuum removed spectral curves of the soil samples.

#### *3.2. Calibration Set and Validation Set Selection Methods*

The MVARC-R-KS method [3] aimed at adaptively selecting the calibration set through multivariate analysis. The procedure of selecting the calibration set using the MVARC-R-KS method comprised two phases [3]. Phase 1 initially mined the clustering distribution zones using the MVARC method, which integrated the a priori algorithm, a density-based clustering algorithm, and the Delaunay triangulation [3]. The association between the SOM content and VNIR spectra may vary according to location and surrounding environment, contributing to the generally strong spatial variability of soil samples. Hence, to open out the spatial heterogeneity of soil samples, environmental variables and SOM content were set as input data for the MVARC method. The spatial heterogeneity of soil samples was checked by Moran's I; results are given in Table 1. Samples in the same cluster mined by the MVARC method were relatively similar with respect to the impact of the surrounding environment. In addition, the influence of environmental variables on SOM in clusters presented significant differences. Phase 2 utilized the Rank–KS algorithm [22] to select the calibration set from the clustering zones. The eventually selected calibration set was an internal component concentration and spectrum representations, as well as an external environment representation. In addition, the size of the calibration set was set as approximately 70% of the soil samples, and the remaining 30% of the samples were used as the validation set, referring to previous studies [23,24].



Apart from the MVARC-R-KS method, three commonly used calibration set and validation set selection methods—-C method, KS method, and Rank–KS method—-were utilized for comparison to verify the representation of the selected calibration set.

#### *3.3. VNIR-Based Prediction Methods*

According to the VNIR-based SOM predicting demands, the PLS–MARS method was proposed to effectively represent the nonlinear relationship between the data with local difference and numerous variables (see Section 3.3.1). To verify the effectiveness of the proposed PLS–MARS method, typical prediction methods, such as MLR, SVM, PLS, and GWR, were used to construct VNIR-based SOM prediction models. To realize a more comparative analysis with the proposed PLS–MARS method, the GWR method was also combined with PLS (named PLS–GWR) to construct a prediction model. Eventually, a series of evaluation indices was used to evaluate the simulation and prediction performance of the prediction models, which is described in detail in Section 3.3.2.

#### 3.3.1. PLS–MARS Method –

The PLS–MARS method is a combination of PLS and MARS methods that constructs a nonlinear prediction model with local regression. The flow of PLS–MARS comprises two parts (Figure 3). Part 1 removes the multicollinearity of the explanatory variables by obtaining PCs using a PLS regression method. Part 2 constructs the MARS model based on the predictors and PCs obtained in Part 1. The details of PLS–MARS are described as follows. – – –

– – **Figure 3.** Procedure of constructing the partial least squares–based multivariate adaptive regression spline (PLS–MARS) model.

<sup>0</sup> <sup>0</sup> PCs were abstracted in Part 1. The explanatory variables and predictors were labeled as *E*<sup>0</sup> and *F*0, respectively. The procedure of obtaining PCs was to abstract the PCs first. Thereafter, a series of PLS regression models based on the obtained PCs and predictors was built. Lastly, the accuracy of the models was evaluated. The model with the highest accuracy was selected, and the corresponding PCs were output as the eventual PCs. The steps of this procedure are described in detail as follows.

<sup>1</sup> <sup>0</sup> <sup>1</sup> <sup>0</sup> Step 1: The PC *t*<sup>1</sup> of the explanatory variables *E*<sup>0</sup> and PC *u*<sup>1</sup> of the predictors *F*<sup>0</sup> were extracted. The estimated values of PCs can be expressed in accordance with Equation (1). The abstracted PCs should represent sufficient information of the variables, and the PCs *t*<sup>1</sup>

and *u*<sup>1</sup> should be highly correlated. The preceding calculation targets imply that *t*<sup>1</sup> and *u*<sup>1</sup> can be calculated based on the conditions in Equation (2) by using the Lagrange multipliers.

$$\begin{cases} \quad t\_1 = E\_0 w\_1\\ \quad u\_1 = F\_0 v\_1 \end{cases} \tag{1}$$

$$\begin{cases} \ \langle t\_1, u\_1 \rangle = \langle E\_0 w\_1, F\_0 v\_1 \rangle = w\_1^T E\_0^T F\_0 \Rightarrow \text{Max} \\ w\_1^T w\_1 = 1, \ w\_1^T v\_1 = 1 \end{cases} \tag{2}$$

Step 2: The regression (i.e., Equation (3)) was constructed based on *t*<sup>1</sup> and *u*1; the model is expressed as follows:

$$\begin{cases} \quad E\_0 = t\_1 \alpha\_1^T + E\_1\\ \quad F\_0 = \mu\_1 \beta\_1^T + F\_1 \end{cases} \tag{3}$$

where *α*<sup>1</sup> and *β*<sup>1</sup> are the model parameters to ensure the high correlation between *t*<sup>1</sup> and *u*1, which can be estimated using least squares criterion; *E*<sup>1</sup> and *F*<sup>1</sup> are the residual matrixes.

Step 3: The residual matrixes *E*<sup>1</sup> and *F*<sup>1</sup> were set as the explanatory variables and predictors, respectively. Step 1 was iteratively repeated until the count of the selected PCs reached *h* (i.e., the matrix rank of *E*0) (Equation (4)). The PCs *t* of the explanatory variables can be calculated using Equation (2).

$$\begin{cases} \quad E\_0 = t\_1 \boldsymbol{\alpha}\_1^T + \cdots + t\_h \boldsymbol{\alpha}\_h^T + E\_h\\ \quad F\_0 = \boldsymbol{\mu}\_1 \boldsymbol{\beta}\_1^T + \cdots \boldsymbol{\mu}\_h \boldsymbol{\beta}\_h^T + F\_h \end{cases} \tag{4}$$

Step 4: The first δ principal component PCs (0 < *δ* ≤ *h*) were successively selected. In addition, the corresponding PLS regression model between the predictors and extracted PCs of the explanatory variables was constructed.

Step 5: The precision of the PLS model was evaluated using the root-mean-square error (RMSEV) calculated by the cross-validation method (leave-one-out sample). High precision of the PLS regression model led to a small value of RMSEV. The PCs *t* = *t*1, · · · , *t<sup>p</sup>* for constructing the PLS regression model with the highest precision were output as PCs for constructing the following the MARS model.

In Part 2, the MARS prediction model was built based on the explanatory variables {*t*1, · · · , *tp*} and predictors. The MARS model was established based on several BFs in the form of an expansion of splines. The estimation of a true function *f*(*t*) utilizing the MARS method [16] based on BFs is expressed as follows:

$$f(t\_1, \dots, t\_p) = \sum\_{m=0}^{M} a\_m B\_m(t\_1, \dots, t\_p) \tag{5}$$

where coefficient *a<sup>m</sup>* is obtained using the least squares method and *B<sup>m</sup> t*1, · · · , *t<sup>p</sup>* (Equation (6)) is one of the BFs, which is multiplied by several *bk*,*<sup>m</sup>* as shown in Equation (7).

$$B\_m(t\_1, \dots, t\_p) = \prod\_{k=1}^{K\_m} b\_{k,m} \left( t\_{v(k,m)} \Big| P\_{k,m} \right) \tag{6}$$

where *K<sup>m</sup>* is the number of *bk*,*<sup>m</sup>* (Equation (7)), which is a bilateral truncation power function decided by the parameters *Pk*,*<sup>m</sup>* and explanatory variables *tv*(*k*, *m*) that correspond to the *k*th truncated spline BF (SBF) in the *m*th term of Equation (6).

$$\left|\boldsymbol{b}\_{\mathrm{k},\mathsf{m}}\left(\boldsymbol{t}\_{\mathrm{v}(\boldsymbol{k},\mathsf{m})}\,\middle|\,\boldsymbol{P}\_{\mathrm{k},\mathsf{m}}\right)\right| = \left[\mathsf{S}\_{\mathrm{k},\mathsf{m}} \times \left(\boldsymbol{t}\_{\mathrm{v}(\boldsymbol{k},\mathsf{m})} - \boldsymbol{r}\_{\mathrm{k},\mathsf{m}}\right)\right]\_{+}^{q} = \max\left(\mathsf{0},\ \mathsf{S}\_{\mathrm{k},\mathsf{m}} \times \left(\boldsymbol{t}\_{\mathrm{v}(\boldsymbol{k},\mathsf{m})} - \boldsymbol{r}\_{\mathrm{k},\mathsf{m}}\right)^{q}\right) \tag{7}$$

where *Pk*,*<sup>m</sup>* = (*Sk*,*m*,*rk*.*m*), truncation direction *Sk*,*<sup>m</sup>* = ±1, *rk*.*<sup>m</sup>* is the knot of the BF *bk*,*<sup>m</sup> tv*(*k*,*m*) *Pk*,*<sup>m</sup>* , and *q* is the power of SBF reflecting the degree of smoothness of the resulting MARS estimation [16]. In this study, *q* was set as 1 to simplify the process.

The process of constructing a MARS model consisted of forward selection and backward pruning. In forward selection, four steps were conducted, as follows:

Step 1: BF was set as *B*0(*t*) = 1, and the threshold of the number of BFs and the threshold of the model precision were set.

Step 2: Two new adaptive BFs that yielded the minimum training error were iteratively created.

Step 3: Step 2 was repeated until the number of BFs reached the number threshold or the model precision reached the threshold.

Step 4: All BFs and related parameters were provided.

In forward selection, permission was limited to adding BFs. Accordingly, several BFs were used merely for constructing succeeding BFs and contributed insignificantly to the eventual model. The threshold of the number of BFs was generally set as high in the forward selection procession. Hence, the backward pruning for deleting redundant BFs was necessary.

In backward pruning, the model was simplified by deleting one least important BF (based on generalized cross-validation (GCV)) at each step until no more BFs were available to be deleted. A new model was rebuilt at each step; the model with the minimum GCV was selected as the eventual prediction model. GCV is computed in accordance with Equation (8).

$$GCV = \frac{\frac{1}{N} \sum\_{i=1}^{N} \left[ y\_i - \hat{f}(t\_i) \right]^2}{\left[ 1 - \frac{M + d \times \frac{(M-1)}{2}}{N} \right]^2} \tag{8}$$

where *N* is the size of the calibration set, *M* is the number of BFs in the model, *d* is a penalizing factor, which is set as 3 in this study according to a report [16,25], and (*M* − 1)/2 denotes the number of the hinge function knots. Consequently, the function penalizes the model for its number of BFs and knots.

In establishing a MARS model, a series of threshold values was preset to realize the adaptive operation and obtain a suitable result. Accordingly, the model with the minimum GCV was selected as the eventual model. The implementation of the PLS–MARS method was realized using MATLAB software.

#### 3.3.2. Fitness Assessment of the VNIR-Based Prediction Model

The prediction model can be evaluated by several assessment indices: the corrected Akaike information criterion (AICc) [26,27], the coefficient of determination of simulation analysis (*R* 2 *c* ), the coefficient of determination of prediction analysis (*R* 2 *p* ), root of mean square simulation error (RMSEV), root of mean square prediction error (RMSEP), relative percent deviation (RPD), and Moran's I. AICc estimates the relative amount of information lost by a given model. When comparing a series of models for the same data, the less information a model loses with a lower AICc value, the higher the quality of that model. *R* 2 *<sup>c</sup>* and RMSEV are used to analyze the simulation precision and stabilization of the model. A low RMSEV value and high *R* 2 *<sup>c</sup>* value indicate the high stabilization and simulation precision of the model. The prediction precision of the model is estimated by *R* 2 *p* , RMSEP, and RPD. If RMSEP is low and *R* 2 *p* is high, then the predictive capability of the model is considered good. In soil spectrographic analysis, if the RRP is less than 1, then the model is not recommended due to poor prediction capability; if the RPD is larger than 1 and less than 1.4, then the prediction capability is still deemed as poor; if the RPD is between 1.4 and 1.8, then the model can be used to perform prediction analysis; if the RPD is larger than 1.8, then the model should have very good prediction capability [28,29]. Moran's I is utilized to test the randomness of residuals. A good model will yield a random series of residuals without autocorrelation [30].

#### **4. Results**

*4.1. Verification of PLS*−*MARS Method Based on a Simulated Data Set*

An illustrative function (Equation (9)) was set to validate the performance of the proposed PLS–MARS method. The method was proposed mainly for nonlinear prediction applications with numerous explanatory variables. Hence, the size of the explanatory variables in the illustrative function in Equation (9) was set as *N* (*N* = 100; *N* = 500), and the corresponding relationship between explanatory variables and the predictor was nonlinear. For comparison, MLR, PLS, and SVM methods were also utilized to estimate the illustrative function. Each prediction model was calibrated by 500 samples generated by the Latin hypercube sampling method [16] and tested using 1000 randomly generated samples.

$$f = \mathbf{x}\_1 \times \sin(\mathbf{x}\_1) + \mathbf{x}\_2 \times \sin(\mathbf{x}\_2) + \dots + \mathbf{x}\_{100} \times \sin(\mathbf{x}\_N) \tag{9}$$

The performance of the models was evaluated by fitness assessment indices (e.g., *R* 2 *c* , RMSEV, *R* 2 *p* , RMSEP, and RPD) (Table 2). The results showed that the MLR model performed well when estimating Equation (9) with 100 explanatory variables, but performed poorly when estimating Equation (9) with 500 explanatory variables. These results suggest the accuracy of the constructed MLR model decreases as the complexity and size of explanatory variables increase. Compared with the MLR model, the PLS method had relatively high accuracy due to the suitability of its equations with numerous explanatory variables. However, SVM and PLS–MARS methods (i.e., as nonlinear prediction methods) had the highest estimation and prediction accuracies when they were used to estimate the relationship in Equation (9).

**Table 2.** Evaluation results of typical prediction models for the illustrative function.


In summary, the comparison experiments on the simulated data set verified that the PLS–MARS method exhibits a high performance for estimating nonlinear relationships to numerous explanatory variables.

#### *4.2. Case Study of PLS–MARS Method*

In this subsection, the PLS–MARS method was used to construct a prediction model between VNIR spectra and SOM content based on the representative calibration set in the study area introduced in Section 2. The study area was utilized as a demonstration zone to validate the effectiveness and accuracies of the PLS–MARS method. Three missions were carried out as follows.

Mission 1 selected the representative calibration set by utilizing the MVARC-R-KS method [3], which is specified in Section 4.2.1. Mission 2 (in Section 4.2.2) constructed the PLS–MARS method, which was calibrated by the selected calibration set using the MVARC-R-KS method [3]. Typical prediction methods based on typical calibration sample selection methods were utilized to predict SOM content in the research area to verify the effectiveness of the proposed PLS–MARS method and the influence of the calibration set. These comparison experiments are elaborated in Section 4.2.3. – –

#### 4.2.1. Calibration Set Selected Using MVARC-R-KS Method

The procedure of selecting the calibration set using the MVARC-R-KS method comprised two phases [3]. Phase 1 consisted of obtaining clustering zones with similar influences of the environment on SOM content. Phase 2 utilized the Rank–KS method to select representative samples from the clustering zones. The selected samples were eventually merged as the calibration set (Figure 4). Typical methods, such as C, KS, and Rank–KS methods, were used for comparison to verify the representation of the calibration set. If the calibration set was representative, then the distribution characteristics of the calibration set would be similar to the entire samples. The statistical values of SOM content and the spectral reflectance information of the calibration sets selected by the aforementioned methods were calculated and are shown in Figure 5 and Table 3. The results showed that the mean and median values of the SOM content in the calibration set selected by the KS method were lower than those in the entire samples. In contrast to the KS method, the C, Rank–KS, and MVARC-R-KS methods obtained calibration sets having a considerably similar distribution to the entire samples for SOM content. Compared with the calibration sets selected by the KS, Rank–KS, and MVARC-R-KS methods, the calibration set selected by the C method had less similar spectral reflectance information to that of the entire samples. Compared with the C, KS, and Rank–KS methods, the MVARC-R-KS method further took into consideration the spatial variation and impacts of the environmental variables [3]. – – – – –

– **Figure 4.** Calibration and validation samples selected utilizing the regional multi-variable associate rule mining and Rank–Kennard-Stone (MVARC-R-KS) method.

– – **Figure 5.** Statistics of the soil organic matter (SOM) content in the following samples: (1) all samples; (2) samples selected utilizing the concentration gradient (C) method; (3) samples selected utilizing the Rank–Kennard-Stone (Rank–KS) method; (4) samples selected using the Rank method; (5) samples selected utilizing the MVARC-R-KS method.

**Table 3.** Statistical values of the spectral reflectance information in the samples.


In summary, the calibration set selected by the MVARC-R-KS method was the SOM content, spectrum representative, and environment representative [3]. Hence, the MVARC-R-KS method seemed to be reasonable to select a representative calibration set.

– 4.2.2. Performance of the PLS–MARS Model Calibrated by the Calibration Set Selected Utilizing the MVARC-R-KS Method

– – – – The PLS–MARS method was used to build the SOM prediction model based on the VNIR spectra in accordance with the selected calibration set (Figure 4). PCs were first calculated and selected. The count of PCs was set at 10, with the highest precision (i.e., RMSEV) of the constructed regression model (Figure 6) based on the PC selection procedure described in Phase 1 of Section 3.2. Thereafter, the PLS–MARS model was constructed based on the PCs and SOM content. The PLS–MARS model was eventually evaluated. Table 4 shows the fit assessments of the prediction model. The result indicated that the PLS– MARS model calibrated by the samples selected utilizing the MVARC-R-KS method [3] could estimate SOM content with high stability and prediction precision.

– –

–

–

–

**–**

–

–

**Figure 6.** Root-mean-square error (RMSEV) of the prediction model constructed using the latent variables (i.e., principal components (PCs)).

**Table 4.** Evaluation results of prediction models calibrated by different calibration set selection methods.


4.2.3. Comparison of the PLS–MARS Model with Typical Prediction Models

The MLR, PLS, SVM, and PLS–GWR methods were utilized to construct the SOM prediction model based on the VNIR spectra in the research region to verify the availability of the PLS–MARS method. The precision of the prediction method was substantially correlated with the selected calibration set. Hence, the C, KS, Rank–KS, and MVARC-R-KS methods were introduced to obtain the representative calibration sets.

The performance of the models was evaluated using fitness assessment indices (e.g., AICc, *R* 2 *c* , RMSEV, *R* 2 *p* , RMSEP, and RPD) (Table 4). Two conclusions could be drawn. One was that the performance of the prediction models highly relied on the selected calibration sets; such performance was achieved for the following reasons. The results imply that the accuracies of the models based on different calibration sets fluctuated heavily. For example, the prediction models calibrated using the calibration sets selected by the C and KS methods had poor prediction capabilities. The Rank–KS method selected the calibration

set for building prediction models with relatively higher accuracies than those of the C and KS methods. The prediction models calibrated by the samples obtained utilizing the MVARC-R-KS method had the highest prediction accuracies. These results sufficiently verified the dependency of the performance of the prediction models on the selected calibration sets. Because the PLS–GWR model was a local prediction model, the calibration sets had less effect on the prediction results. Another conclusion was that the PLS–MARS method outperformed the typical prediction models for estimating VNIR-spectrum–based SOM content in the study area. The following results support this conclusion. The results in Table 4 show that the MLR model had poor prediction accuracies and thus could not be used to predict SOM content in the study area. The possible reason is due to the inapplicability of these methods for complex high-dimension data sets. Although the PLS and SVM models had relatively higher simulation and prediction accuracies, their model accuracies were lower than the PLS–GWR and PLS–MARS models that considered the local difference. The Moran's I indexes in Table 1 indicate that both the PLS–MARS model and PLS–GWR model could maintain the spatial heterogeneity distribution characteristics of data. In addition, random spatial residuals without autocorrelation could further verify the effectiveness of the PLS–GWR model and the PLS–MARS model. Compared with the PLS–GWR model, the PLS–MARS model, which considered the nonlinear relationship between SOM content and VNIR spectra, had the highest simulation and prediction accuracies among the prediction models. Since PLS–GWR and PLS–MARS were constructed based on the same PCs, their AICc values were comparable. Nevertheless, the AICc values indicated that the proposed PLS–MARS model had higher model quality than the PLS–GWR model. In summary, the PLS–MARS model calibrated by the samples obtained utilizing the MVARC-R-KS method indicated immense potential for building a prediction model of SOM content based on the VNIR spectra in the riverside region.

#### **5. Discussion**

Intrinsic strategies for the prediction method and selected calibration set were significant factors that affected the performance of the constructed prediction model. The influence of the calibration set on model performance, and the good performance of the MVARC-R-KS method and the proposed PLS–MARS prediction method will be discussed in Sections 5.1 and 5.2, respectively.

#### *5.1. Influence of the Calibration Set on the Model Performance*

If the calibration set can well represent the relation of the SOM content with VNIR spectra, then the prediction model calibrated using this calibration set will have a good prediction performance. Table 4 shows that the prediction models calibrated using the calibration sets selected by the classical C and KS methods lacked good prediction performance, whereas the prediction models based on the Rank–KS method could build models with considerably high *R* 2 *<sup>p</sup>* and RPD in the study area. This phenomenon was mainly due to the fact that the calibration set obtained utilizing the C or KS method was representative of either the SOM content or the spectrographic information. By contrast, the calibration set obtained utilizing the Rank–KS method was both SOM content and spectral information representative. These results also indicated that the component concentrations and spectrographic information should be simultaneously taken into account to obtain the representative calibration set to accurately build a VNIR-spectrum–based SOM prediction model. However, in comparison with the Rank–KS method, the MVARC-R-KS method could select the more representative calibration sets, which could calibrate the prediction models with the highest accuracy. This result could be attributed to the consideration of the impact of the surrounding environment and spatial heterogeneity on the samples in the MVARC-R-KS method, which were substantially consistent with real situations.

In summary, the prediction models depended substantially on the selected calibration sets. A representative calibration set instead of a randomly selected calibration set should be obtained in advance to construct a prediction model with good performance. Hence,

calibration set selection should be taken seriously during the construction of the prediction models. This study suggests that the MVARC-R-KS method [3] can be set as the optimal choice in the riverside region. This is because this method simultaneously considers the influence of multiple influential factors (i.e., chemical components, spectrographic information, and environmental factors), thereby making the VNIR-spectrum–based SOM prediction model calibrated by the selected calibration set an extensively applicable one.

#### *5.2. Strategies of the PLS–MARS Method*

The performance of the PLS–MARS method was analyzed using the simulated data set and real application. In addition, the effectiveness and accuracy of PLS–MARS were compared with those of several commonly used prediction methods (i.e., MLR, PLS, SVM and PLS–GWR methods) in literature. The results in the simulated application implied that the PLS, SVM, and PLS–MARS methods could efficiently estimate the relationship between numerous explanatory variables and their responses. However, if the relationship exhibited significantly non-linear features and spatial heterogeneity, such as the relation of the SOM content with VNIR spectra in the study area, then the PLS–MARS method produced more accurate results than the rest of the prediction methods tested. This superiority is generally due to the consideration of the natural local difference of data sets, and the construction of a local non-linear optimization model in the PLS–MARS method rather than the global non-linear optimization characteristics of the other prediction methods. The adaptive strategy of the PLS–MARS method also facilitated the automatic estimation of natural relationships and supported their promotion. Thus, the PLS–MARS method has immense potential for estimating the complex nonlinear relationships of geographical data with numerous explanatory variables, and the PLS–MARS method calibrated by the samples obtained utilizing the MVARC-R-KS method shows great potential in estimating significantly nonlinear relationships in spatial and non-spatial domains.

#### **6. Conclusions**

The PLS–MARS method was proposed to construct a nonlinear model for SOM content prediction, in which the MVARC-R-KS method was adopted to select the optimal calibration set. Several commonly used calibration set selection methods were applied to calibrate the PLS–MARS model, and their influences on the prediction precision of the established model were investigated. The implementation procedure for the proposed method was described in detail. The proposed approach was illustrated and validated through a simulated data set and practical application in the Jianghan Plain. A comparative study with conventional prediction methods indicated that the proposed PLS–MARS method presents the following advantages: (1) the proposed method can efficiently estimate the nonlinear relationship underlying data with numerous variables; (2) the proposed method can accurately construct the highly nonlinear relationship of the SOM content with VNIR spectra, with consideration of the spatial heterogeneity of the relationship; (3) the proposed method can automatically construct a prediction model without substantial prior knowledge of the data set.

The following new findings were obtained based on the application of the PLS–MARS method to estimate SOM content in the Jianghan Plain in China: (1) the accuracy of prediction models, including the PLS–MARS model, depends significantly on the selected calibration sets; (2) the combined PLS–MARS and MVARC-R-KS methods have immense stability and prediction capabilities when estimating the relation of SOM content with VNIR spectra in the study area.

In general, the proposed PLS–MARS method can accurately estimate the nonlinear relationship with local differences. The PLS–MARS method calibrated by the samples obtained utilizing the MVARC-R-KS method has considerable potential for estimating SOM. Furthermore, the PLS–MARS method is conceptually simple and easily executed via a programming language, thereby ensuring easy application.

Future studies will focus on further practical applications of the PLS–MARS method. For example, a PLS–MARS model calibrated by the calibration set selected utilizing

MVARC-R-KS method can serve as a potential prediction model to estimate other soil properties (e.g., iron content, copper content, and soil organic C) based on the VNIR spectra in the study area or other areas. The proposed PLS–MARS method can also be used to characterize the nonlinear relationships of other geographic phenomena.

**Author Contributions:** X.W. and M.Z. conceived and designed the experiments; X.W. and C.Y. performed the experiments; all the authors analyzed the data; X.W. wrote the paper; all authors contributed to the revision of the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by [the Scientific Research Program of the Hunan Education Department] grant number [19C1113], [the Natural Science Foundation of Hunan Province] grant number [2020JJ5363], and the National Natural Science Foundation of China grant number [41901314].

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data available on request due to restrictions privacy. The data presented in this study are available on request from the corresponding author. The data are not publicly available because the data is the result of the whole research group's hard work.

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

#### **References**


#### *Article*

## **Orchard Floor Management A**ff**ects Tree Functionality, Productivity and Water Consumption of a Late Ripening Peach Orchard under Semi-Arid Conditions**

**Pasquale Losciale 1,\* , Liliana Gaeta <sup>2</sup> , Luigi Manfrini <sup>3</sup> , Luigi Tarricone <sup>4</sup> and Pasquale Campi <sup>2</sup>**


Received: 23 October 2020; Accepted: 13 November 2020; Published: 17 November 2020 -

**Abstract:** Semi-arid conditions are favorable for the cultivation of late ripening peach cultivars; however, seasonal water scarcity and reduction in soil biological fertility, heightened by improper soil management, are jeopardizing this important sector. In the present two-year study, four soil managements were compared on a late ripening peach orchard: (i) completely tilled (control); (ii) mulched with reusable reflective plastic film; (iii) mulching with a Leguminosae cover-crop flattened after peach fruit set; (iv) completely tilled, supplying the water volumes of the plastic mulched treatment, supposed to be lower than the control. Comparison was performed for soil features, water use, tree functionality, fruit growth, fruit quality, yield and water productivity. Even receiving about 50% of the regular irrigation, reusable reflective mulching reduced water loss and soil carbon over mineralization, not affecting (sometimes increasing) net carbon assimilation, yield, and fruit size and increasing water productivity. The flattening technique should be refined in the last part of the season as in hot and dry areas with clay soils and low organic matter, soil cracking increased water evaporation predisposing the orchard at water stress. The development and implementation of appropriate soil management strategies could be pivotal for making peach production economically and environmentally sustainable.

**Keywords:** mulching; flattening; irrigation; photosynthesis; transpiration; soil quality; water stress integral; fruit growth; water use efficiency; productivity

#### **1. Introduction**

Fruit growing is a key sector for the Mediterranean economy, society and environment. It is the highest value among the agricultural productions, representing 17% of the total EU agricultural turnover (FAO—Food and Agriculture Organization of the United Nations, 2018). Furthermore, orchards contribute to land preservation via stewardship and climate regulation through evapotranspiration [1] and represent one of the most typical fruit crops of the Mediterranean Basin [2], thus being at the basis of both its economy and dietary culture.

Fruit growing has a double-faced nature. On one side, the great demand of high-quality products. The fruits, most of which are delivered to the fresh market, are asked to meet the consumer demand with very high-quality standards [3] as fruit consumption improves health and well-being [4]. This, at

world level, makes the fruit sector very competitive, creating concerns and often pushing the fruit growers more on the yield than in quality production. Orchard intensification techniques are one of the results from this request. On the other side, society has expressed concern about the exploitation of agricultural inputs because of their dramatic impact on natural resources and ecosystem functioning [5].

Peach is among the most representative and valued fruit species in the Mediterranean Basin. The southern Italy environment is usually hot and dry until the beginning of autumn, therefore is particularly suitable for early and late ripening peach cultivars. However, late cultivars are more water demanding due to the long-lasting persistence of fruits on the plant [6].

Climate change and foreseeable future resource limitations (mainly water) threaten Mediterranean fruit production. The secure supply of high-quality fruit is under jeopardy because of increased heat and water stress conditions [7], reduced productive land [8], decreased soil fertility due to intensive management practices, reduced water availability and competition for water with other productive sectors and human activities [8]. In Mediterranean countries, the reduction in organic matter, which impacts soil fertility, is a common occurrence, for climatic reasons [9]. This is often aggravated by improper management of the residues of agricultural products, the low use of organic amendments and the rapid mineralization of the organic compounds due also to intensive tillage practices [10]. A multi-year life cycle assessment study showed that fertilizers and energy consumption (i.e., electricity, fossil fuel) and water were the main impact factors in peach cultivation [11]. Fertilizers and energy consumption are indeed indicated to be the most impactful on emissions and climatic disorders. Climate change, with its increased temperatures, can have negative effects on tree productivity if scheduling irrigation is not applied properly. However, agriculture is already the main user of water, consuming about 70% of freshwater, and it must achieve savings, rather than further increases, in its water needs [12]. Water scarcity can be tackled by improving water saving techniques. In particular the rational management of soil can increase water use efficiency either through decreasing soil evaporation, using artificial [13–15] or natural [16,17] mulching material, and increasing soil water holding capacity, via decompaction and organic matter enrichment. Recent experiences in horticultural crops [18,19] reported that water and N recycling at agroecosystem level can be enhanced by cover crop practice and natural mulching covering techniques, independently from the soil management strategy. The authors pointed out that the improved nitrogen surplus was not sufficiently retained in the agroecosystem without cover crop. On the contrary, the practices adopted in the treatments with the cover crop or temporary intercropping considerably improved the N self-sufficiency of the system. Beneficial effects have been found in grapes [20,21]. Grape ecosystem services provided by Mediterranean vineyards are particularly threatened, because soil functions are often impaired by yearly repeated intensive agricultural practices or weed and pest management. The authors demonstrated that the potential of soil management practices to enhance soil functioning, can be promoted by the presence of a cover crop, even temporarily, in the inter row [20,21]. Moreover, Almagro found that improved soil management in rainfed Mediterranean agroecosystems can be a powerful strategy to mitigate the current atmospheric CO<sup>2</sup> increase, through soil carbon sequestration and stabilization [22]. However, few publications on the effect of orchard floor modification on fruit trees, in general, and on peach tree development and functioning, including fruit production, in particular, are available [23]. These authors found that changed soil management practices such as zero tillage, supply of organic amendments, understory mowing, retention of crop residues, can result in worthwhile gains within a long-term period. In a Mediterranean peach orchard, these gains can include increased production and also increased sustainability with higher level of soil organic carbon and litter carbon pools [23]. However, to date, orchard soil and water management often rely only on grower and extension service experience or, in the most advanced cases, are driven by data on soil water content and/or climate conditions and plant status [24–26]. Real-time tree performance, as well as the inter-relation among the different chemical, physical and microbiological variables affecting soil fertility, is little considered. Solutions to the major threats encountered by the fruit production sector may be found through approaches that consider the soil–plant system as a whole and, therefore, address improving the entire orchard performance to

cope with water scarcity and climate change [27,28]. Since orchards are particularly complex systems, different methods and approaches should be adopted under different, even contrasting, pedo-climatic, economic and social conditions [29]. Even lesser explored in horticulture, the use of different mulching material in the fruit orchard is deserving of particular interest. The use of plastic mulch was adopted with the aim of increasing the water use efficiency in a dryland rainfed area. However, the plastic adopted in the mulching was dark in color with negligible results on light diffusion and probably enhancing soil temperature [30,31]. Recent studies pointed out the positive effect of mulching with high-reflective biodegradable plastic film on productivity and water use efficiency on peach [15].

The aim of the present study was to investigate the effect of four different orchard floor managements on tree functionality, productivity and water consumption of a late ripening peach cultivar orchard under semi-arid conditions. The comparison was performed among the following treatments: (i) completely tilled, (ii) mulched with reusable and reflective plastic film, (iii) mulching with a Leguminosae cover crop flattened after peach fruit set, with the aim of increasing the organic matter and water holding soil capacity; (iv) completely tilled and reducing irrigation at the same volumes supplied to the plastic mulched treatment.

#### **2. Materials and Methods**

#### *2.1. Experimental Set-Up and Pedo-Climatic Conditions*

The trial was carried out in 2015 and 2016 at the Experimental farm of the Council for Agricultural Research and Economics (Research Centre for Agriculture and Environment), in southern Italy (Rutigliano, lat.: 40◦59′ N, long.: 17◦01′ E, alt.: 147 m asl) on 3-year-old peach (*Prunus persica* (L.) Batsch var. *laevis*) trees of a late ripening cultivar "Calred" [32] grafted on Missour rootstock, trained as slender spindle and spaced 4.0 × 2.5 m. The experimental site is under the Mediterranean climate, characterized by warm and dry summers. The average air temperature throughout the year and during the vegetative–reproductive season is 15.5 and 20 ◦C, respectively, and the annual rainfall is about 535 mm, mainly concentrated in the autumn and late winter periods and usually greatly reduced, or absent, in the spring–summer period [33].

Four different orchard managements were tested: soil tilled (T); inter-row mulching with a reusable and machine-resistant reflective plastic film (C/820 Black Silver Orchard; thickness: 100 µm; Ginegar Plastic product Ltd., Ginegar, Israel) to reduce soil evaporation and to increase the diffuse light (M); inter-row mulching with horse bean (*Vicia faba.* L.) sown in November and flattened after peach fruit set forming a natural mulching on the soil (F). The last treatment (S) was established on tilled soil supplying the same irrigation volume of M that was supposed to be lower than the control (T) as evaporation was limited by the plastic mulching. Since horse bean contributed to Nitrogen fixation, N supply on F treatment was halved in comparison to T, M and S, while the Phosphorus and Potassium supply was increased by 25% in order to feed the service crop.

In order to verify the homogeneity of soil characteristics at the beginning of the trial, as well as to evaluate the evolution of the soil conditions as a function of the four treatments, three soil samples per treatment were collected at the beginning of the trial (before flattening period), after the harvest of the first year and at the end of the second year (just before the winter pruning). Soil was evaluated for its physico-chemical traits: soil texture by hydrometer method, total carbon organic content (TOC, %) by the dry-combustion procedure with a TOC Vario Select analyzer (Elementar, Germany), pH, electric conductibility (EC, dS m−<sup>1</sup> ), N (g kg−<sup>1</sup> , Kjeldahl procedure) and P (mg kg−<sup>1</sup> , Olsen method) content. Soil texture was similar among the four treatment and it was classified as clay loam [34]. Soil water content in volume at field capacity (FC, −0.03 MPa) and wilting point (WP, −1.5 MPa) were 0.34 m<sup>3</sup> m−<sup>3</sup> and 0.21 m<sup>3</sup> m−<sup>3</sup> , respectively (measured using the Richards chambers). At 0.6 m of depth, the parent rock is present; this reduces the capacity of the root systems to expand beyond this layer. At the beginning of the trial also, the evaluated chemical trait results were not statistically different among the treatments (Table 1).


**Table 1.** Soil chemical traits at the beginning of the trial (14 April 2015) after the harvest of the first year of the study (24 September 2015) and before the winter pruning of the second year (16 December 2016). Within each date and for each variable, different letters indicate a statistical difference at *p* ≤ 0.05.

#### *2.2. Water Supply and Soil Water Content*

Water was supplied with drip irrigation system having 2 drippers per tree, and a flow rate of 8 l h−<sup>1</sup> per dripper. Volumetric soil water content (SWC) was measured by capacitive probes (10HS, Decagon Devices Inc., Pullman, WA, USA) linked to dataloggers (Grillobee, TecnoEL, Italy). For each treatment, three points were monitored. Capacitive probes were installed horizontally into the soil profile, on the row at 0.3 m from peach trees and at −0.1, −0.3 and −0.5 m from the soil surface, in order to intercept the dynamics of soil water content below the dripping lines. For each treatment, water content in the whole soil profile was calculated averaging the values of the three depts in each of the three points. Probes were previously calibrated in order to measure the volumetric soil water content (SWC) and identify the intervention threshold (IT). The IT corresponded to the SWC at which the readily available water was completely used. The IT of 0.26 m<sup>3</sup> m−<sup>3</sup> was adopted; this value was obtained considering a depletion fraction (fraction of available soil water that can be depleted from the root zone before moisture stress) of 0.5 [35]. When the IT was reached, the amount of water necessary to return at FC was supplied [36]. T, M and F were irrigated monitoring the SWC while S received the same water volume of M.

#### *2.3. Leaf Functionality and Tree Water Relations*

Three plants similar in canopy size and potential crop load were selected for each treatment. At fruit cell division, pit hardening, fruit cell expansion and close to the harvest stages, leaf net photosynthesis (Pn, µmol m−<sup>2</sup> s −1 ), stomatal conductance (gs, mol m−<sup>2</sup> s −1 ) and transpiration (Tr, mol m−<sup>2</sup> s −1 ) were measured on well-exposed leaves placed on the east and west side of the canopy, 4 times during the day (9.00–17.00 h), with an open circuit infrared gas analyzer fitted with an LED light source (Li-COR 6400XT, LI-COR, Lincoln, Nebraska, USA). At each time of the day and canopy side, light intensity was maintained constant across the 4 treatments, setting the LED light source at the natural irradiance experienced by the leaf immediately before the measurement. The values obtained on the west and east side of the canopy were averaged for each tree. At the same time of measure, stem water potential (Ψs, MPa) was measured on the same trees belonging to the 4 treatments according to [37]. Pn and Tr, collected during the day, were integrated [38], providing the specific amount of CO<sup>2</sup> (ΣPn, mol m−<sup>2</sup> ) and water (ΣTr, mol m−<sup>2</sup> ) fixed and transpired by a square meter of leaf during the time of measure. To take into account the amount of energy used by trees to raise water from the soil during the time of measure, the water stress integral (SΨ, GPa) was calculated as the difference between the integral by time of stem water potential and the integral by time of the minimum stem water potential measured in the same time range [39,40].

#### *2.4. Fruit Growth and Productivity*

The fruit growth pattern was monitored during the season by means of a digital caliper implemented with a datalogger able to store the data (HK—Horticultural Knowledge s.r.l. Bologna, Italy) on twelve fruit per tree. The fruit volume (V, cm<sup>3</sup> ) and the absolute growth rate (AGR, cm<sup>3</sup> day−<sup>1</sup> ) were calculated assuming the shape of the peach as a spheroid and measuring the three axes of each peach [15].

At harvest (ready for hand selection: when fruit flesh firmness was around 3–5 kg cm−<sup>2</sup> ), the number of fruits per tree (NF), the average fruit weight (FW, g), the yield (Y, t ha−<sup>1</sup> ) and the irrigation water productivity (WPi, kg) of fresh fruit per cubic meter supplied with irrigation [41] were evaluated on the same trees monitored for fruit growth and leaf functionality. The total soluble solids content (TSS, ◦Brix), flesh firmness (FF, kg cm−<sup>2</sup> ) and the percentage of fruit skin red overcolor (RC, %) were measured on 10 fruit per tree.

#### *2.5. Statystical Analysis*

For each period of measure soil data, ΣPn and ΣTr were subjected to ANOVA. Fruit growth data (V and AGR) used for the statistical analysis were obtained averaging the measures taken on each tree. A by-time repeated ANOVA was performed analyzing separately the data of fruit cell division, pit hardening, fruit cell expansion and ripening stages, respectively. Three productivity, water use efficiency and fruit quality variables were tested by means of an ANCOVA considering the number of fruits as the covariate variable.

#### **3. Results**

#### *3.1. Pedo-Climatic Conditions*

The two-year study's thermic patterns during the experiment period (1 June–10 September) were almost comparable with an average minimum, maximum and average temperature of about 20.6, 29.5 and 25.0 ◦C, respectively. The year 2015 was less rainy than 2016 with a cumulative rainfall during the period of 126 mm. The year 2016 showed a cumulative rainfall of 206 mm till 31 August and an additional 147 mm of rain fallen in the first fortnight of September (Figure 1).

Soil measurements were performed on 14 April 2015, 24 September 2015 and 16 December 2016, before peach's full bloom, after the harvest and before the winter pruning of the second season, respectively. At the beginning of the trial, the soil chemical traits evaluated resulted in being not statistically different among the four treatments (Table 1). N, P, pH, conductibility and total organic carbon content were about 1.0 g kg−<sup>1</sup> , 37.2 mg kg−<sup>1</sup> , 8.4, 0.16 dS m−<sup>1</sup> and 1.0 g kg−<sup>1</sup> , respectively. After the first harvest, the treatments continued to be similar for soil chemical traits. At the end of the two-year trial, S showed the lowest value of P and TOC while the highest P content was recorded for F, followed by T and M (Table 1).

#### *3.2. Water Supply and Soil Water Content*

In order to have all the treatments at the same soil moisture conditions at the beginning of the experiment, three and two full irrigations were provided regardless of the treatments in 18 May, 29 May and 5 June in 2015, as well as in 6 June and 14 June, in 2016. In 2015, the seasonal water supply was 1557 m<sup>3</sup> ha−<sup>1</sup> for T and F, and 815 m<sup>3</sup> ha−<sup>1</sup> for M, S. The irrigation season lasted about

4 months (Figure 2A). In 2016, T and F received 1810 m<sup>3</sup> ha−<sup>1</sup> while M and S 1023 m<sup>3</sup> ha−<sup>1</sup> and the duration of the irrigation season was about 3 months (Figure 2B). M and S received about 50% less water than T and F in both years. Soil water content remained between the field capacity (FC) and the intervention threshold (IT) till the end of July in T, M and F during the two seasons. In the same period, S showed SWC lower than IT for several days (Figure 3, June–July). From the beginning of the irrigation season till the end of July, M and F showed the highest SWC values (average SWC of 0.32 and 0.33 m<sup>3</sup> m−<sup>3</sup> in 2015 and 2016, respectively, for M; 0.31 and 0.32 m<sup>3</sup> m−<sup>3</sup> for F). In the same period, T revealed SWC similar to M and F, in 2015 (average of 0.31 m<sup>3</sup> m−<sup>3</sup> ), while in 2016, it decreased at an average value of 0.29 m<sup>3</sup> m−<sup>3</sup> . The average SWC of S in the same period was 0.29 and 0.26 m<sup>3</sup> m−<sup>3</sup> in 2015 and 2016, respectively (Figure 3). Due to local watershed restrictions occurring during the period August–September for both the years, the irrigation frequency (number of peaks in Figure 3) was reduced for all the treatments. In this period, the SWC of T, F and S fell below the intervention threshold several times and in S it reached the wilting point; M showed the highest soil water content, rarely below IT (Figure 3). In this period, the comparison between T and F, receiving almost the same water volume in each irrigation, showed that in August–September, after water supply SWC declined faster in F than in T, reaching values lower than T and close to WP at the end of the irrigation season (Figure 3). − − − −

**Figure 1.** Air temperature (lines) and rainfall (bars) recorded for the peach orchard under investigation during the two years of study.

− − **Figure 2.** Cumulated water volumes supplied to T, F (closed circles) and to M, S (open triangles) recorded in 2015 (**A**) and 2016 (**B**).

−

#### − − *3.3. Leaf Functionality and Tree Water Relations*

Full bloom occurred on 7 April and 14 March in 2015 and 2016, respectively. Leaf gas exchange and stem water potential measures were performed 55, 86, 87, 88 and 120 days after full bloom (DAFB) in 2015 (1/6, 2/7, 3/7, 4/7 and 5/8), and 92, 107, 120, 135, and 163 DAFB (14/6, 29/6 12/7, 27/7 and 24/8) in 2016.

**A B** 

− **Figure 3.** Soil water content (m<sup>3</sup> m−<sup>3</sup> ) pattern recorded on the 4 treatments in 2015 (**A**) and 2016 (**B**). Horizontal dotted and continuous lines represent the field capacity (FC), the wilting point (WP) and the intervention threshold (IT), respectively.

3.3.1. Season 2015

Σ Σ <sup>Ψ</sup> − At 55 DAFB, when water supply differentiation was not established yet, cumulative photoassimilation (ΣPn), transpiration (ΣTr) and the integral water stress (SΨ) were similar among the four treatments (Figure 4; Table S1). The average air temperature (Tair) and vapor pressure deficit (VPD) during the measure period (9:00–17:00) were 29.9 ◦C and 2.3 kPa, respectively; the midday stem water potential was around −0.9 MPa and the average stomatal conductance was about 0.164 mol m−<sup>2</sup> s −1 . At 86–87 and 88 DAFB, trees were at the end of the pit hardening stage and water supply was differentiated among the four treatments. At 86 DAFB, the average VPD and air temperature were 3.5 kPa and 36.4 ◦C, respectively. A slight reduction in ΣPn and ΣTr were observed in S, while the water stress integral SΨ was statistically lower in S and F in comparison with T and M (Figure 4; Table S1). The midday stem water potential was around −1.2 MPa in M and T and −1.5 MPa in F and S. The average gs recorded during the measure period was about 0.102 mol m−<sup>2</sup> s −1 for M and T, 0.094 mol m−<sup>2</sup> s −1 for F and 0.085 mol m−<sup>2</sup> s −1 for S. At 87 DAFB, the average VPD was about 4.5 kPa and the air temperature recorded during the time of measure was about 37.5 ◦C. S showed the cumulative photoassimilation and transpiration to be lower than the remaining treatments (0.142 and 72.5 mol m−<sup>2</sup> , respectively). The highest S<sup>Ψ</sup> was recorded in M (−3.6 GPa), followed by T (−7.0 GPa); F and S showed the lowest

values of the water stress integral of about −10 GPa (Figure 4; Table S1). The midday stem water potential was −1.1, −1.25, −1.3 and −1.4 MPa in M, T, F and S, respectively, while the average gs was around 0.105 mol m−<sup>2</sup> s −1 in M, T, F and 0.064 mol m−<sup>2</sup> s −1 in S. At 88 DAFB, the average VPD and the air temperature of the period of measure (9:00–17:00) were 3.8 kPa and 37.1 ◦C, respectively. The highest cumulative photoassimilation was recorded on M (0.25 mol m−<sup>2</sup> ), followed by T and F (0.198 and 0.188 mol m−<sup>2</sup> , respectively). S revealed the lowest ΣPn (0.125 mol m−<sup>2</sup> ) and the same trend was observed for the cumulative transpiration (Figure 4A,B; Table S1). The average gs recorded during the period was 0.134 mol m−<sup>2</sup> s −1 for M followed by T and F, with an average gs of about 0.096 mol m−<sup>2</sup> s −1 and S (gs, 0.057 mol m−<sup>2</sup> s −1 ). M had the highest S<sup>Ψ</sup> (−2.5 GPa) followed by T (−6.7 GPa); the lowest S<sup>Ψ</sup> values were recorded in F and S with values of −14.2 and −17.2 GPa, respectively (Figure 4C). At 120 DAFB, trees were in the full fruit cell expansion stage. The average (9:00–17:00) air temperature and VPD were 33.5 ◦C and 2.5 kPa, respectively. The lowest ΣPn and ΣTr were recorded in S (0.166 and 48.6 mol m−<sup>2</sup> , respectively) while the remaining treatments were similar (Figure 4A,B; Table S1). The same trend was observed for SΨ reaching the lowest levels of the season (Figure 4C; Table S1). The midday stem water potential was −1.4, −1.5, −1.6 and −1.8 MPa for M, T, F and S, respectively, and the average stomatal conductance recorded within the measure period (9:00–17:00) was around 0.106 mol m−<sup>2</sup> s −1 for M, T, F, and 0.066 mol m−<sup>2</sup> s −1 for S.

#### 3.3.2. Season 2016

At 92 DAFB, when water supply was the same for all the treatments, no differences were recorded in terms of cumulative photoassimilation, transpiration and water stress integral; the average VPD and air temperature during the period of measure (9:00–17:00) were about 2.0 kPa and 29.4 ◦C, respectively; the midday stem water potential was −0.6 MPa and the average gs was about 0.182 mol m−<sup>2</sup> s −1 . At pit hardening (107 DAFB), ΣPn, ΣTr and S<sup>Ψ</sup> were similar among the treatments (Figure 5; Table S1). Average air temperature and VPD were 29.6 ◦C and 2.4 kPa, respectively; the midday stem water potential was −0.8 MPa and the average gs was 0.115 mol m−<sup>2</sup> s −1 . At the beginning of fruit cell expansion (12 July, 124 DAFB), the average VPD was about 4.5 kPa and air temperature 37.6 ◦C. T and F showed the highest ΣPn (around 0.23 mol m−<sup>2</sup> ), while the lowest one was recorded on S (∼0.15 mol m−<sup>2</sup> ). M revealed an intermediate ΣPn of 0.18 mol m−<sup>2</sup> (Figure 5A; Table S1). The highest cumulative transpiration was recorded on F (146.1 mol m−<sup>2</sup> ), followed by T (127.0 mol m−<sup>2</sup> ); the lowest ΣTr was observed in M and S with an average cumulative transpiration of 89.8 mol m−<sup>2</sup> (Figure 5B). The average gs during the period of measure was about 0.107 mol m−<sup>2</sup> s −1 for T and F and 0.073 mol m−<sup>2</sup> s −1 for M and S. F and T had a quite similar water stress integral (−3.9 GPa), higher than S<sup>Ψ</sup> recorded in S and M of about −12 GPa (Figure 5C). The midday stem water potential was −1.0 MPa for T and F, and −1.3 MPa for S and M. At the fruit cell expansion stage (135 DAFB), the average air temperature and vapor pressure deficit recorded during the period of measure were about 32.5 ◦C and 2.8 kPa, respectively. M and F showed ΣPn and ΣTr higher than T and S (Figure 5A,B; Table S1). The average gs recorded from 9:00 to 17:00 was about 0.118 mol m−<sup>2</sup> s −1 in M and F, and 0.098 mol m−<sup>2</sup> s −1 in T and S. The water stress integral was −7.8 GPa in T, followed by M (−10.6 GPa); S and R had the lowest S<sup>Ψ</sup> of about −14.8 GPa (Figure 5C; Table S1). The midday stem water potential at 135 DAFB was about −1.0 MPa in T, −1.1 MPa in M and −1.3 MPa in F and S. Close to the harvest (163 DAFB), the average Tair and VPD were 29.1 ◦C and 1.9 kPa, respectively. M showed a cumulative net photoassimilation and transpiration higher than T and S and the lowest values were observed in F (Figure 5A,B; Table S1). The average stomatal conductance followed the same trend with gs of 0.126 mol m−<sup>2</sup> s −1 in M, about 0.099 mol m−<sup>2</sup> s −1 in T and S, and 0.064 mol m−<sup>2</sup> s −1 in F. S<sup>Ψ</sup> in M was −4.4 GPa, higher than T and S (~−9.0 GPa); the lowest water stress integral was observed in F with S<sup>Ψ</sup> of −13.4 GPa (Figure 5C). The midday stem water potential followed the same trend with values of −1.4, −1.6, −1.6 and −1.8 MPa recorded in M, T, S and F, respectively.

Σ Σ Ψ **Figure 4.** Cumulative leaf net photosynthesis (ΣPn) (**A**), transpiration (ΣTr) (**B**) and water stress integral (SΨ) (**C**) calculated for T (black bars), M (grey bars), F (dashed bars) and S (dotted bars) during the time of measure (9:00–17:00 h) of each day of measurement in 2015. Within the same date different letters indicate a statistical difference at *p* ≤ 0.05.

Σ Σ Ψ **Figure 5.** Cumulative leaf net photosynthesis (ΣPn) (**A**), transpiration (ΣTr) (**B**) and water stress integral (SΨ) (**C**) calculated for T (black bars), M (grey bars), F (dashed bars) and S (dotted bars) during the time of measure (9:00–17:00 h) of each day of measurement in 2016. Within the same date different letters indicate a statistical difference at *p* ≤ 0.05.

≤

#### *3.4. Fruit Growth and Productivity*

#### 3.4.1. Season 2015

At the end of the fruit cell division stage (38–77 DAFB), no differences for fruit volume were recorded among the four treatments while the average AGR values observed in M and S were lower than those measured on T and F (Table 2; Figure 6). During the pit hardening stage (77–105 DAFB), the fruit volume was similar among the treatments (Table 2), and the absolute growth rate was higher in T, M and F (~0.59 cm<sup>3</sup> day−<sup>1</sup> ) than in S with an AGR value of 0.43 cm<sup>3</sup> day−<sup>1</sup> (Table 2). The reduced AGR in S was observed starting from 94 DAFB (Figure 6B). In the full fruit cell expansion stage

− −

(105–125 DAFB), T, M and F continued to have fruits bigger than S (Table 2) with a difference between S and the remaining treatments growing progressively (Figure 6A). The average AGR recorded between 105 and 125 DAFB was higher in M, T and F (~1.07 cm<sup>3</sup> day−<sup>1</sup> ) than in S with an AGR of 0.73 cm<sup>3</sup> day−<sup>1</sup> (Table 2). F maintained an AGR similar to M and T till 115 DAFB; afterwards it decreased, reaching values closer to S (Figure 6B). During the last days before the harvest (125–148 DAFB), M and T showed an average fruit volume (~105.8 cm<sup>3</sup> ) higher than F and S with a value of about 88 cm<sup>3</sup> (Table 2). The same behavior was observed for the absolute growth rate with values of about 2.34 cm<sup>3</sup> day−<sup>1</sup> , for M and T and 1.86 cm<sup>3</sup> day−<sup>1</sup> for S and F (Table 2).

**Figure 6.** Fruit growth (**A**) and absolute growth rate (**B**) pattern recorded in 2015 for T, M, F and S treatments.

**− − −**

≤


**Table 2.** Average fruit volume (V) and absolute growth rate (AGR) recorded for T, M, F and S treatments in different fruit growth stages (Days After Full Bloom range, DAFB Range), in 2015. Within each stage, a by-time repeated ANOVA was performed. For each variable, different letters indicate a statistical difference at *p* ≤ 0.05.

In 2015, the third leaf season for the peach orchard, the harvest was performed in two picks: on 3 and 7 September, 149 and 153 DAFB, respectively. A late fruit drop occurred during the season modifying the initial crop load imposed. The ANCOVA revealed the significative effect of the number of fruit as the covariate variable for the yield (F = 63.26; *p* < 0.001), fruit weight (F = 19.98; *p* = 0.001) and water productivity (F = 30.98; *p* < 0.001), whose values have been adjusted accordingly. The yield was similar among the treatments (~11.00 t ha−<sup>1</sup> ), while the fruit fresh weight was higher in M and T (166.9 and 149.4 g, respectively) than in F and S with values of 112.9 and 93.6 g, respectively (Table 3). Water productivity was higher in M and S (average of 12.01 kg m−<sup>3</sup> ) than in T and F with an average value of 6.93 kg m−<sup>3</sup> (Table 3). No differences for the total soluble solid content was observed, recording an average value of about 18 ◦Brix. The percentage of fruit skin over color was higher in M (~92.3%) than in the remaining treatments (~56.7%) while the flesh firmness was higher in F, followed by T, S and M with values of 3.7, 3.0, 2.8, and 2.2 kg cm−<sup>2</sup> , respectively (Table 3).


**Table 3.** Yield (Y), fruit fresh weight (FW), water productivity (WPi), sugar content (TSS) flesh firmness (FF) and fruit skin red overcolor, measured on the 4 treatments under investigation in 2015. Data were subjected to ANCOVA analysis, considering the number of fruits per tree as a covariate variable and were adjusted accordingly. For each variable, different letters indicate a statistical difference at *p* ≤ 0.05.

#### 3.4.2. Season 2016

During the pit hardening (57–107 DAFB) and the first part of fruit cell expansion stages (107–140 DAFB), no differences for fruit volume and AGR were recorded among the four treatments (Table 4; Figure 7). Afterwards (140–164 DAFB), the average fruit volume for the period was similar among the treatments, while the absolute growth rate was higher in M and T (5.89 and 5.16 cm<sup>3</sup> day−<sup>1</sup> , respectively) than in S and F, with values of 4.63 and 4.31 cm<sup>3</sup> day−<sup>1</sup> , respectively (Table 4). Starting from 150 DAFB, a divergent pattern for AGR was observed comparing T and M versus F and S: while in the former treatments AGR continued to increase, in the latter ones it decreased (Figure 7B). Close to the harvest, the fruit volume of M and T was higher than F and S (Figure 7A).


**Table 4.** Average fruit volume (V) and absolute growth rate (AGR) recorded on T, M, F and S treatments in different fruit growth stages, in 2016. Within each stage, a by-time repeated ANOVA was performed. For each variable, different letters indicate a statistical difference at *p* ≤ 0.05.

**Figure 7.** Fruit growth (**A**) and absolute growth rate (**B**) pattern recorded in 2016 on T, M, F and S treatments.

≤ **− − −** In 2016, the harvest occurred on 26 August (165 DAFB) with a single pick. Even in this season a late fruit drop affected the imposed crop load. The ANCOVA showed a significant effect of the number of fruits as a covariate variable for yield (F = 126.95; *p* < 0.001) and WPi (F = 76.84; *p* < 0.001), whose values have been adjusted accordingly. Y was similar for the four treatments, with values ranging from 16.78 t ha−<sup>1</sup> (M) to 13.56 t ha−<sup>1</sup> (F). Fruit fresh weight was higher in M (256.65 g) than in T (234.46 g) and the lowest values were observed in S and F with FW of 216.92 and 211.65 g, respectively (Table 5). The highest WPi was recorded in M (16.2 kg m−<sup>3</sup> ), followed by S (14.03 kg m−<sup>3</sup> ), T (8.38 kg m−<sup>3</sup> ) and F with WPi of 7.45 kg m−<sup>3</sup> (Table 5). No difference in terms of sugar content, flesh firmness and percentage of red overcolor was observed among the treatments (Table 5).


**Table 5.** Yield (Y), fruit fresh weight (FW), water productivity (WPi), sugar content (TSS) flesh firmness (FF) and fruit skin red overcolor, measured on the 4 treatments under investigation in 2016. Data were subjected to ANCOVA analysis, considering the number of fruits per tree as a covariate variable and were adjusted accordingly. For each variable, different letters indicate a statistical difference at *p* ≤ 0.05.

#### **4. Discussion**

The thermic pattern during the two years of study was quite similar and in line with the average of the site while year 2015 was less rainy than 2016 (Figure 1). The similar chemical–physical conditions of the soil recorded at the beginning of the trial suggested that the starting point for the four treatments was the same (Table 1). At the end of 2016, the different floor and water managements affected the soil chemical features (Table 1). S, subjected to water restriction and tillage, showed a decrease in TOC and assimilable P. As observed in other studies, under semi-arid conditions the high temperature and evapo-transpirative demand of the environment, jointly with the low soil moisture lasting for long time during the year, may have increased the mineralization processes, reducing the already low organic carbon in the soil [10,42,43]. An additional, related consequence could have been the reduction in P in its assimilable form in the tilled treatment subjected to water shortage [44]. A decrease in P in its assimilable form was observed in an evergreen forest, in the Mediterranean area, subjected to drought stress [45]. M, experiencing the same water supply of S, did not show a decrease in TOC suggesting the positive role of artificial, reflective mulching in preventing organic matter over degradation (Table 1). The use of organic mulching (F) did not affect TOC, while it was not possible to verify if the increase in P occurring in this treatment was given by the floor management per se or to the increased amount of P supplied as fertilizer to feed the service crop (Table 1). Artificial reflective mulching reduced soil evaporation behaving as a physical barrier against the water loss but also reducing solar radiation absorption by the soil and the wind speed at the surface [46]. Even when receiving about 50% less water than the control (T), M maintained SWC higher than T for most of the irrigation season in both years (Figure 3). These preliminary results confirmed what was observed in other studies on rain-fed peach orchards where mulching with plastic film reduced water loss of about 15% in comparison with tilled soil [14]. The flattening technique (F) ameliorated SWC conditions in the rainier year (2016) till the end of July. Afterward, its SWC dropped faster than T after irrigation, in both the years (Figure 3). This rapid water depletion might have been attributable to the disruption of the natural mulching and to the onset of relevant soil cracking occurring in clay soils, low in organic matter [47]. August was characterized by a general water limitation caused by watershed restriction usually occurring in this area (water demand for crop and for civil use increases in this period). Under this stressful condition, M continued to maintain SWC higher than the remaining treatments (Figure 3); only in a few cases its soil water content dropped below the IT, assuring an adequate soil moisture for all the season long [13,15]. When the soil evaporation was not contrasted, the reduction in water supply (S), produced a progressive consumption of the readily available water with values of SWC very low and below the IT (Figure 3).

Soil and water management strategies affected leaf functionality and water relations in peach trees. In both the years, when water supply was not yet differentiated among the four treatments and SWC was within the readily available water range (Figure 3), T, M, F and S behaved similarly (Figures 4 and 5) for carbon assimilation (ΣPn), water transpiration (ΣTr) and plant water status, expressed as water

stress integral (SΨ). Reflective mulching, preventing the excessive soil evaporation and increasing the diffuse light, maintained the pedo-climatic conditions favorable to photosynthetic activity of leaves [13,14] for the entire vegetative–reproductive season in the two years of study (Figures 4A and 5A). Excluding the measure performed in 2016 at 124 DAFB (12 July), ΣPn, ΣTr and S<sup>Ψ</sup> in M were similar, and in some cases higher than T receiving double the amount of water (Figures 4 and 5). When water shortage was associated to tillage (S), tree water status (SΨ) was affected and net photosynthesis was subjected to stomatal limitation (reduction in average gs and ΣTr) and probably to non-stomatal limitation as well (Figures 4 and 5). Previous research in pear, apple, peach and grapevine suggested that stomatal closure could affect leaf thermoregulation inducing the increase in leaf temperature and raising the activity of photorespiration [48–50]. This process is a photoprotective strategy for the plant, but it means a loss of carbon in terms of biomass accumulation [51–53]. The use of natural mulching associated with full irrigation (F) did not affect leaf functionality (ΣPn and ΣTr) in comparison with T till the end of July, in both years (Figures 4 and 5). This suggested the absence of such competition between the main and the service crop and the positive effect of the flattening technique in controlling weeds. This positive effect was already observed in vegetable crops [54,55]. In August, when the flattened crop was disrupted and the soil cracking caused a rapid water depletion (Figure 5, 163 DAFB), SΨ, ΣPn and ΣTr values in F were the lowest of the four treatments. From a leaf functioning point of view the use of the flattening technique for all of the dry season appeared to be detrimental under particular pedoclimatic conditions such as hot summer, clay and poor of organic matter soils.

The differences among the treatments appeared more evident when high temperature and VPD were associated with soil moisture limitation. The comparison between the days 86 and 87 DAFB of year 2015 revealed that SWC did not change markedly within each treatment as well as SΨ (Figure 4C). However, at 87 DAFB, the environment was more water demanding than the previous day: the average Tair and VPD passed from 36.4 to 37.5 ◦C and from 3.5 to 4.5 kPa, respectively. As a consequence, soil with an adequate SWC allowed leaves to maintain the stomata opened, increasing ΣTr and ΣPn. On the other hand, S, having a low SWC, reduced gs and leaf transpiration [56], thus carbon assimilation (Figure 4A). Trees, being in the middle of the Soil Plant Air Continuum (S.P.A.C.) were strongly influenced by the status of rhizosphere and air. The comparison between F and S at 86, 87 and 88 DAFB revealed that, even if the two treatments had the same SΨ and midday stem water potential, S showed a cumulative net photosynthesis and leaf transpiration lower than F (Figure 4). The same behavior was observed in 2016, comparing M and S at 124 DAFB and F and S at 135 DAFB. This late ripening peach cultivar seemed to have a "pessimistic" (also called conservative or near iso-hydric) behavior. At low soil water content and high vapor pressure deficit, S sustained its stem water potential, at the same level of F, above a safe threshold to prevent embolism [57,58]. This defense strategy was at the expense of CO<sup>2</sup> fixation since it was regulated by stomatal closure [59,60]. These findings are in contrast with those described by Xiloyannis et al. (1980) on another late ripening peach cultivar considered aniso-hydric [61], suggesting the needing to deepen this issue.

During the first year of study, the pattern of fruit growth of the four treatments was similar till the end of fruit cell division stage (38–77 DAFB); however, M and S, receiving less water, showed an average absolute growth rate lower than T and F (Table 2). This difference was completely recovered in M during the pit hardening stage (77–105 DAFB), while S continued to have an average AGR lower than the remaining treatments (Table 2). Passing from pit hardening to the fruit cell expansion stage (105–125 DAFB), peach fruit became more water and carbon demanding [62]. Water shortage, jointly with the reduction in net photosynthesis, led to fruit growth limitation. As revealed by the average V and AGR in the period 105–125 DAFB (Table 2; Figure 6), the reduced AGR initially observed in S during the pit hardening (77–105 DAFB) resulted in fruit size lower than the remaining treatments (Table 2). In August, during the last part of fruit cell expansion and close to the harvest (125–148 DAFB), fruit growth was limited even in F, suggesting that the rapid water depletion occurring in this period affected fruit volume and growth rate (Table 2, Figure 6). All the advantage of F on S was lost, and, at the end of the season, the two treatments showed the smallest fruits (Figure 6, Table 3). The

same behavior was observed in 2016. The rainier season alleviated the effect of water shortage and till 140 DAFB (first part of fruit cell expansion) no differences were recorded among the treatments for fruit volume and AGR (Table 4, Figure 7). During the last part of fruit cell expansion (140–164 DAFB), even in this season, F and S showed the lowest average AGR (Table 4), fruit volume (Figure 6) and fruit size (Table 5).

Fruit yield was generally low considering that Calred is a late ripening cultivar. However, it should be taken into account that this cultivar was subjected to a late fruit drop in both the years and that, probably, trees were not at a fully mature productive status, as in 2015 and 2016, they were at third and fourth leaf, respectively. The different orchard floor management and water supply did not affect yield, but the fruit size with biggest fruit picked on T and M (Tables 3 and 5). Water productivity of M was about 70 and 90% higher than T suggesting that the use of artificial reflective mulching could be considered a water friendly strategy in rainfed [14] as well as in irrigated peach orchards (Tables 3 and 5). Although in 2015 peaches seemed to ripe earlier in M in comparison with the remaining treatments (Table 3), fruit quality (TSS, FF and RC) was generally not affected by the different managements (Tables 3 and 4).

#### **5. Conclusions**

Under semi-arid conditions and where water supply is limited, the choice of an appropriate orchard floor management could be of pivotal importance for getting the peach production both economic and environment-friendly. Even receiving about 50% of the regular irrigation, reusable reflective mulching reduced water loss and soil carbon over mineralization, not affecting (sometimes increasing) net carbon assimilation, yield, and fruit size. As a consequence, water productivity was drastically increased. These first results suggested that the reflective mulching strategy could be considered to be water and soil "friendly". This management technique is firstly described and explored on a peach orchard, thus the studies on the development and the use of alternative material for mulching should be explored. The flattening technique as mulching strategy should be refined for the final part of the irrigation season, especially in those hot and dry areas with clay soils, low in organic matter, thus predisposed to cracking.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2076-3417/10/22/8135/s1, Table S1: Anova results for Cumulative leaf net photosynthesis (ΣPn), transpiration (ΣTr) and water stress integral (SΨ) calculated in each day of measurement in 2015 and 2016. For each variable, within the same date the asterisk indicates *p* ≤ 0.05.

**Author Contributions:** Conceptualization, P.L.; data curation, L.G., L.M., L.T. and P.C.; formal analysis, L.G. and P.C.; investigation, L.G., L.T. and P.C.; methodology, P.L. and P.C.; supervision, P.L.; validation, P.L. and L.G.; writing—original draft, P.L. and L.M.; Writing—review and editing, P.L., L.G. and L.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to thank Di Gennaro D. and Amendolagine A.M. for their contribution in data collecting, and Introna P. and Volpicella M. for their valuable operative effort in conducting the orchard.

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

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Review* **Bioremediation of PAH-Contaminated Soils: Process Enhancement through Composting**/**Compost**

**Tahseen Sayara <sup>1</sup> and Antoni Sánchez 2,\***


Received: 22 April 2020; Accepted: 22 May 2020; Published: 26 May 2020

**Abstract:** Bioremediation of contaminated soils has gained increasing interest in recent years as a low-cost and environmentally friendly technology to clean soils polluted with anthropogenic contaminants. However, some organic pollutants in soil have a low biodegradability or are not bioavailable, which hampers the use of bioremediation for their removal. This is the case of polycyclic aromatic hydrocarbons (PAHs), which normally are stable and hydrophobic chemical structures. In this review, several approaches for the decontamination of PAH-polluted soil are presented and discussed in detail. The use of compost as biostimulation- and bioaugmentation-coupled technologies are described in detail, and some parameters, such as the stability of compost, deserve special attention to obtain better results. Composting as an ex situ technology, with the use of some specific products like surfactants, is also discussed. In summary, the use of compost and composting are promising technologies (in all the approaches presented) for the bioremediation of PAH-contaminated soils.

**Keywords:** bioremediation; composting; PAHs; organic co-substrates; soil

#### **1. Introduction**

Globally, different anthropogenic activities have resulted in increasing environmental pollution, and its consequences has injured almost all components of the ecosystem [1–3]. Soil, as a vital component of the terrestrial ecosystem, is prone to pollution from different sources, including industrial and agricultural activities [4–7]. Wide verities of pollutants entering the soil posing a huge threat and risk to human health and natural ecosystem [8–12]. Polycyclic aromatic hydrocarbons (PAHs), petroleum, and related derivatives represent the main sources of soil contamination [13–17]. Indeed, these organic pollutant groups are listed as priorities and receive considerable attention, owing to their toxic, genotoxic, mutagenic, and potentially cancer-causing properties [18,19].

To deal with this problem, several treatment technologies are used, including chemical, physical, and biological, as well as thermal for remediation of these contaminated soils. Among the best approaches is the bioremediation technology, which is categorized as a promising approach that continues to gain more attention due to its efficiency, cost-effectiveness, and environmental-friendly byproducts [20–22]. The process mainly relies on the activity of a wide spectrum of microorganisms to degrade the target contaminants to lower toxic levels. Bioremediation of PAH-contaminated soil has been performed utilizing distinctive approaches [7]. In any case, composting as a remediation approach has been considered a reasonable strategy in this field, because it provides nutrients for indigenous microorganisms to degrade the target contaminants; simultaneously, applying this approach is a great opportunity for feasible and sustainable reuse of the natural biodegradable fraction of wastes. Additionally, the process is cost-effective compared with other approaches—for instance, composting

costs between \$50–\$140 per ton, while applying slurry or biopiling treatments cost \$170 per ton and \$130–\$260 per cubic meter, respectively [23–27]. Bioremediation of PAH-contaminated soil through composting could be implemented through incorporating PAH-contaminated soils to the composting process, or by adding compost to contaminated soils. Also, bioaugmentation or surfactant application might be included to achieve the final set objectives [16,25,28–34]. Biodegradation of PAHs intrinsically depends on microbial activity, where bacteria and fungi are considered the foremost vital variables governing the bioremediation process [35–37]. However, the functionality of these microorganisms is affected by different factors within the composting mixture, including biotic and abiotic factors. In this context, the environmental condition (pH, temperature, moisture,), nutrient availability, oxygen presence, and bioavailability of the contaminants are essential parameters for process control and performance [38].

This review focuses on the application of composting and compost addition for the bioremediation of soils contaminated with PAHs. In this regard, the impact of different controlling factors like temperature, PAH structure and concentration, co-substrate stability, co-substrate mixing ration, and bioaugmentation are discussed. Moreover, other issues, such as bioavailability, surfactant application, and the degradation pathways of PAHs are illustrated, in order to provide an insight into the process that is necessary for new development.

#### **2. Soil Contamination with PAHs**

Soil represents a vital component of all terrestrial ecosystems. However, it is subjected to degradation or decline in its quality as a result of different anthropogenic activities that have resulted in increasing the rate of contamination [4,5,7]. Therefore, polluted land is a worldwide concern, and can be viewed as major obstruction to sustainable development and modern environmental protection [39]. Soil contamination has been recognized as one of the major dangers to soil function in Europe by the Communication from the European Commission "Towards a Thematic Strategy for soil Protection" [40,41]. The issue has expanded with expanding public awareness and concern about the presence of chemicals in the environment, particularly due to their different unfavorable impacts on the ecosystem and human health. Polycyclic aromatic hydrocarbons (PAHs) have been recorded as pollutants of priority importance due to their properties and ubiquitous occurrence, as well as their recalcitrance [18,42,43]. Consequently, great efforts worldwide have been directed toward remediating these pollutants from the environment.

#### *2.1. PAHs: Properties and Sources*

PAHs are a group of ubiquitous organic pollutants with at least two aromatic rings (Figure 1), and are poorly soluble in water (Table 1). Due to their chemical structure, PAHs have hydrophobic properties, which refers to their ability to accumulate on the surface of solid materials like soil, sediment, sewage sludge, and solid wastes. The dangers emerging from the presence of PAHs in soil are related to the toxic nature of those pollutants [42,43]. It is noteworthy that some substances in this group have been recognized as mutagenic, carcinogenic, and teratogenic [18,19].

Sources of PAHs are categorized into natural as well as anthropogenic sources: hydrothermal process volcanoes, forest fires, and waste burning are natural sources of PAHs. Anthropogenic sources include waste incinerators, burning of fossil fuels during heating processes, incomplete combustion of organic matter, petrochemical spills on land, wood burning, petrol and diesel oil combustion, gasification, and plastic waste incineration [13–15,17]. Globally, 16 to 32 PAH compounds are subjected to mandatory control, due to their harmful properties [44]. PAH persistence and hydrophobicity in environmental components are the main factors that exacerbate the pollution problem, taking into account that soils receives a considerable share of this pollution (sink), due to their complex matrix structure that facilitates the sorption of these pollutants. Soil organic matter is a decisive factor in determining the degree of PAH sorption into the soil, along with the physicochemical properties of

the PAHs themselves [18,19,23,45,46]. Therefore, the remediation of soils polluted by aged PAHs has become a major issue for environmental scientists in recent years [12–14,47,48].

**Figure 1.** Structure and chemical formula of the 16 polycyclic aromatic hydrocarbons (PAHs) listed as priority pollutants by the United States Environmental Protection Agency (USEPA).

#### *2.2. Bioremediation of PAH-Contaminated Soils*

When natural biodegradation processes cannot achieve the desired goals, in this case, human intervention becomes necessary to stimulate the process above naturally occurring microbial process [49]. Accordingly, several approaches have been used to enhance bioremediation efficiency. These approaches, which could be used separately or in combination (two or more) include, but are not limited to, biostimulation (providing nutrients for increasing the microbial activity), bioaugmentation (introducing a consortium of indigenous or exogenous microorganisms), using surfactants, and co-metabolism [50,51]. Recently, various studies have been done in attempt to understand the

process hierarchy and to provide solutions for different process limitations. For instance, much research has been carried out to better understand the microbial behavior and its interaction with the contaminants during the bioremediation process, whereas others have focused on introducing exogenic and genetically engineered microbes for process enhancement [52].


**Table 1.** Selected properties of the 16 USEPA PAHs.

\* **Kow:** octanol–water partition.

#### **3. Composting Technology**

Composting is defined as an aerobic process, which fundamentally requires oxygen, optimal moisture content, and porosity to stabilize the organic waste, and its common control variables are temperature, oxygen, and moisture [53]. Thus, composting bioremediation is the adaptation and application of the composting technology for wastes and contaminant treatments.

In order to achieve optimum results within a reasonable time during any composting treatment, process-controlling parameters have to be adjusted within the optimum values, and the process passes through two main stages. First is the decomposition/active stage, which is characterized by extensive microbial activity that leads to a steadily increase in the temperature, passing from the mesophilic ranges (25–45 ◦C) to reach the thermophilic ones (more than 45 ◦C). To maintain aerobic conditions for effective microbial activity during this stage, a high rate of aeration is needed. Second is the curing stage; this take place at a lower temperature, and microbial activity is relatively low, as the nutrients pool has been depleted. Material humification is an important characteristic occurring in this stage [54], which gives an interesting value to the produced compost, especially for soil bioremediation, as will be discussed later in this work.

#### **4. Bioremediation of PAH-Contaminated Soil by Composting**

Composting technology is categorized as ex situ technology, which has been used for the treatment of contaminated soils. During the last few decades, the process received more attention, as it has proved its high efficiency in degrading various organic contaminants like, among others, PAHs, pesticides, explosives, and chlorophenols [25,55–60]. Essentially, the process relies on the addition of compost or organic co-substrates/amendments to the contaminated soil, and while the co-substrate matures, due to the action of various microbial populations within the mixture, the target pollutants

are degraded [57,61]. Thus, treatment of PAH-contaminated soil combined with composting of organic waste could be an interesting option and a sustainable method with much increasing attention. It would enable eco-friendly disposal of such waste and enhance the biodegradation rate of PAHs [7,60,62]. The biodegradation process efficiency depends fundamentally on the bioavailability of the substrates, environmental conditions (pH, moisture, temperature), the presence of oxygen, and the availability of nutrients [38]. Remarkably, the bioremediation of PAH-contaminated soils through composting has confirmed this technique's capability to overcome most obstacles that might hinder reaching its goal, which is the removal of contaminants [10,63–67].

As the process is based on mixing the contaminated soil with organic co-substrates, any failure may result in producing much greater quantity of contaminated material, and this is recognized as the main concern of using such approach. This weakness, and the general scarcity of information on the toxicity, distribution, and bioavailability of such contaminants in compost-amended soils, may therefore result in the drawing up of excessively stringent soil assessment measures with remediation cost implications [68].

#### *4.1. E*ff*ect of PAH Characteristics and Concentrations*

The physical and chemical properties of PAHs have a considerable effect on their biodegradation rate. Microbial assimilation and biodegradation of these compounds basically depends on their solubility. Nevertheless, most of compounds belonging to this group are characterized as poorly soluble in water, especially with their increasing molecular weight and angularity (Table 1, Figure 1), which thus increase their hydrophobicity [69–72]. This was obvious in many studies dealing with the biodegradation of different PAHs. For instance, Han et al. [73] investigated the application of different agricultural waste on the biodegradation of aged PAHs in soil microcosms over 90 days. The initial concentration of total PAHs in the soil was 36.1 mg kg−<sup>1</sup> dry soil, where four-ring PAHs comprised 41.7% of the total PAHs. The results demonstrated higher degradation rates of 40.7–61.2% for PAHs with low molecular weight (LMW), compared to 18.7–33.1% for those with high molecular weight (HMW) in all soil microcosms. Similarly, Luki´c et al. [74] showed that LMW-PAH removal was more favorable in the mesophilic phase, with 11% and 15% residues in the soil, than in the thermophilic phase, with 29% and 31% residues. Additionally, more resistance to degradation was observed for HMW PAHs, resulting in a decrease in the total removal, which was less than 50% for both benzo[a]pyrene and benzo[k]flouranthene, in all treatments [75,76]. In this regard, even though both compounds have the same number of benzene rings (five) and molecular weights, the higher octanol–water partition coefficient (log *Kow*) of benzo[k]flouranthene increased its hydrophobic properties and consequently its degradation rate under the same conditions. Indeed, higher log *Kow* leads to a higher potential of bioaccumulation, which is the main factor responsible for the lower biodegradability of such compounds [77]. Obviously, and according to the obtained result in different studies, there is a consistent relationship between the persistence of PAHs in the environment and increasing their numbers of benzene rings, which ultimately affects their biodegradation rate.

The concentrations of PAHs also have a substantial influence on the microbial activity in such treatments, since high concentrations would lead to toxic or inhibition conditions. Meanwhile, low concentrations could be below the rate needed to stimulate microbial cultures to degrade these contaminants [78,79]. This was obvious in the study conducted by Sayara et al. [78], in which the PAH concentrations had a crucial effect. Low concentrations were found to be less than the rates that are assumed to initiate the degradation process, since microbial communities prefer the utilization of readily available nutrients, which are consumed quickly before initiating biodegradation of the target PAHs. The same results were obtained by Zappi et al. [80], where low concentrations of PAH did not degrade, even when the system was supplanted with additional carbon sources. Wu et al. [66] showed that compost addition is an effective approach for enhancing PAH removal from soils, but increasing the ratio of added compost does not necessarily help to increase removal. Nevertheless, enhanced removal by compost addition seems more effective for higher initial PAH concentrations. In

this regard, Jorgensen et al. [81] demonstrated that the degradation rate of a compound is proportional to its concentration, especially for highly soluble compounds, and argued that the degradation of hydrocarbons is governed by first-order kinetics. However, this argument may be validated to some extent, as high concentration may become detrimental to microbial activity and disturb nutrient balance, especially when LMW PAHs are present [78].

#### *4.2. E*ff*ect of Temperature*

Providing optimum temperature is an intrinsic factor for the successful biodegradation of PAH. The importance of this factor stems from its influence on the metabolic activity, bioavailability, solubility, and diffusion rate of the target contaminate [82]. It is noteworthy that the solubility of PAHs increases with temperature, which ultimately increases the bioavailability of the PAH molecules. However, increasing temperature is associated with decreasing oxygen solubility, which on turn reduces the metabolic activity of aerobic microorganisms. Furthermore, and to a certain extent, the specified temperature range will determine the types of dominant microorganisms and their enzymatic activity that will undertake the degradation [73].

The successive stages during the normal composting process (mesophilic phase, thermophilic phase and curing phase) are expected to be accompanied by specific populations of bacteria, and different effects on contaminants are found with different stages of compost product. The biodegradation of PAHs occurs over a wide temperature range, and microorganisms have found to be adapted to biodegrade PAHs at extreme temperature conditions. Under mesophilic and thermophilic temperatures ranges, it has been found that the enzymatic activity of microorganisms increases, which helps in increasing the rate of hydrocarbon degradation. However, it should be underlined that a great amount of research has been directed to focus on the process behavior under mesophilic conditions, as it is believed that a wide spectrum of microbial communities could be present in active roles under these temperatures, and thus reasonable degradation rates could be achieved [78,83–86].

During in-vessel composting of pyrene-contaminated soil, composting temperature affected the prevailing of some microbial groups over others, and the predominant bacterial community changed over time. The degradation of pyrene was dominated by α-, β-, and γ-Proteobacteria, as well as *Actinobacteria*, at 38 ◦C during 14 days of composting, and then *Streptomyces* at 55 ◦C. Later, at 70 ◦C and after 42 days of composting, *Acinetobacter* and *Thermobifida* occupied leading position. Finally, *Thermobifida* and *Streptomyces* flourished after 60 days of composting at 38 ◦C [87]. Concerning the temperature effect, Luki´c et al. [74] claimed that degradation rates of 89% and 59% for three-ring and four-ring PAHs, respectively, were achieved in reactors under mesophilic temperatures. In contrast, reactors displaying a thermophilic range ended with 71% and 41% removal for the same pollutants, respectively, during the bioremediation process. The addition of compost significantly promoted the removal of PAHs and alkanes up to 88% after 50 days of incubation under mesophilic temperatures (28 ◦C), compared to the natural biodegradation of hydrocarbons in soils without compost [67]. Additionally, the composting of PAH-contaminated soils under different conditions and different organic substrates were found to perform better under mesophilic conditions [23,25,78,88,89]. LMW-PAH concentrations, such as naphthalene, acenaphthylene, acenaphthene, fluorene, anthracene, and phenanthrene, were decreased by an average of 89% at 38 ◦C, which is twice that compared to the concentration reduction at 55 ◦C, which was an average of 45%. Simultaneously, no big difference was observed concerning HMW PAHs, including fluoranthene, pyrene, benzo[a]anthracene, and chrysene, where the removal rate was by an average of 67% at 38 ◦C, compared to 69% at 55 ◦C. Nevertheless, a high temperature was considered adverse to microbial activity, and volatilization was the leading mechanism of PAH removal [88]. Under these conditions, it is assumed that a longer incubation period under the mesophilic phase could facilitate PAH removal, due to the richest microbial diversity and possible increased microbial activity [27]. According to these studies, and others in the literature, mesophilic temperatures demonstrated their viability and were found to be more favorable for degrading LMW PAHs, with great success in many cases due to the large microbial diversity; however, these temperatures were not found to be so efficient in the degradation of recalcitrant PAHs [69,75,88,90].

On the other hand, thermophilic ranges have been documented as enhancing PAH degradation. During the composting of hydrocarbon-polluted sediments (total petroleum hydrocarbons (TPH) = 40.3 gkg−<sup>1</sup> dw) with different organic co-substrates, Alves et al. [91] point out that fish sludge achieved higher temperatures and was able to maintain thermophilic temperature longer than other amendments, which ultimately led to greater TPH removal rates (39.5%). It was assumed that such conditions are conducive to develop fungal communities and exert a surfactant effect, thus promoting the degradation rates. Similarly, Zhu et al. [92] proposed that enhanced solubility under thermophilic conditions could explain the higher removal rate (46%) of benzo[a]pyrene in composting treatments compared to 29% under mesophilic ones. However, whether the increased solubility or microbial community changes contribute to the high-temperature impacts needs further investigation. Viamajala et al. [82] further demonstrated that the elevated temperature during the thermophilic phase of composting enhanced the solubilization rates of phenanthrene, and hence its degradation. Based on the aforementioned observations, it is clear that the impact of composting temperature is correlated to the physiochemical properties of the targeted PAH, as the corresponding degrading microorganism are specific to temperature [93]. Generally, and despite of the different observations, mesophilic temperatures and the dominant microorganisms under these conditions are believed to be more preferable for the degradation of such compounds [69,78,83–86,94].

#### *4.3. E*ff*ect of Organic Co-Substrate Stability*

Even though various organic co-substrates/amendments have demonstrated their viability in the composting of PAH-contaminated soils, composition of these materials varies significantly in the sources and stages of decomposition [25,59,78,95,96], which as a consequence influences the removal rate in different ways [59,69,78,97]. The selected organic co-substrates for the bioremediation process should contribute in improving and overcoming any deficiencies or limitations that influence the process performance and efficiency. Accordingly, selection of the most suitable organic co-substrates represents as a major challenge in such studies [59,95].

In the bioremediation of PAH-contaminated soils, organic matter stability is of particular importance, as this parameter is directly correlated to the organic substrates' composition and biological activity [98]. Various studies have pointed out that the fate of PAHs is dependent on the quality and nature of the amended organic matter [25,69,78,93]. Bioremediation of PAH-contaminated soils with more stable compost has proved to be more effective than with less stable or fresh organic amendments [25,78,97]. In this context, the preference of these substrates related to the presence of humic substances was found to form a considerable part of stable compost and was proportional with its degree of stability [25,78]. This humic matter was documented to enhance the organic compounds' bioavailability [78,97]. During the composting process, organic co-substrates provided nutrients for microorganisms [99]; meanwhile, humic matter evolution is expected to facilitate the microbial accessibility to PAHs. This behavior is established as a result of decreasing humic matter binding affinity and increasing of the heterogeneity of binding sites, closer to soil humic matter, which is conducive to microbial accessibility to PAHs [97,100]. Additionally, stable compost contains low biodegradable organic matter content and a higher concentration of organic macromolecules [101], which are believed to enhance the biodegrading of the contaminant. The presence of easily degradable organic matter is assumed to reduce the process efficiency, as microbial cultures prefer to use easily degradable organic matter and thus decrease or retard utilization of the contaminant. Another important point in this item is represented by effect of the potential working surface area. In this regard, less degraded organic matter generally has coarse fractions (>5 mm), whereas most humified organic matter is generally presented in fine fraction [102]. The finest compost size fraction (<3 mm) with a higher surface area ratio provides more accessibility to microorganisms and releases more nutrients compared to coarse compost fraction [101,103].

During the bioremediation process, an increase in the content of humic matter from 0.23% to 0.70% was observed, and these changes resulted from the structural changes that occurred in the material composition. Potentiometric titrations of humic acid solution showed increases in the buffering and redox capacities of humic acids [104]. Plaza et al. [97] reported that the composting process caused a structural conversion of humic acids from an organic substrate by reducing the aliphatic fraction and increasing the polarity and aromatic polycondensation in a PAH-contaminated soil. This conversion decreased the PAH binding affinity of humic acids, and thus improved PAH-degrading microbial accessibility. Similar results were observed in other studies, which supports the application of stable organic co-substrates [24,66,73,78,105,106].

It is worth mentioning that when the same ratio of inorganic fertilizer (Nitrogen (N) .Phosphorus (P) .Potassium (K)) was compared with organic compost on the bioremediation of diesel-polluted agricultural soil over a two-month period, the results revealed that total petroleum hydrocarbon removal from polluted soil was 71.40 ± 5.60% and 93.31 ± 3.60% for N.P.K. and compost-amended options, respectively [107]. Also, after 30 weeks, the removal efficiencies of TPH in the soils were 29.3% under natural attenuation, 82.1% when nutrients (NH4NO<sup>3</sup> and K2HPO4) were added, and 63.7% when the mixture was supplemented with 20% (*w*/*w*, dry weight basis) of aged refuse. However, a removal efficiency of 90.2% was recorded when nutrient and aged refuse were combined together. Nutrients plus aged refuse made the TPH concentration decrease to below the threshold level of commercial use required for Chinese soil quality for TPH (<3000 mg/kg) in 30 weeks. It was also found that dehydrogenase activity, bacterial counts, and degrader abundance in the soil were remarkably enhanced by the addition of aged refuse (20% *w*/*w*) [108]. All these results confirm the suitability of stable compost over other organic and inorganic substrates. Therefore, one can conclude that introducing an adequate organic co-substrate is usually more efficient in enhancing the bioremediation process, as observed in different studies. This advantage presumably corresponds to the capacity of the compost to perform simultaneously for both bioaugmentation and biostimulation.

#### *4.4. E*ff*ect of the Mixing Ratio*

The suitability of different substrates based on their physiochemical properties is recognized as an important factor in the composting of PAH-contaminated soils. Determining the appropriate quantity to be added to the mixture is also of great importance, since an inappropriate ratio may hamper or inhibit microbial activity and bioavailability [78,93]. It has been determined that even though microbial metabolism may be temporarily increased using a certain mixing ratio, the long-term inhibition of functionally important organisms may result in the failure of the bioremediation of high-molecular-weight PAHs [78]. The amount of various nutrients, and the ratio of nutrients like C, N, and P in particular, are quite conceivable as being involved in the success of the bioremediation process. Furthermore, determining the minimum quantity of the amendment that could support and maintain the desired activity with a high degradation rate is directly related to process economics [78].

As reported in the study conducted by Wang et al. [109], a microorganism's selection of nutrients could delay the degradation of pollutants, as normally microorganisms prefer easily degradable materials over resistant ones. This study revealed that that treatments with amendment ratios of 1:1 and 2:1 had average TPH removal rates of 30.7% and 33.3%, respectively, but the amendment ratio of 3:1 had a slower net degradation rate of between 11.6% and 26.8%. An excess of readily degradable carbon might overtake the TPH and act as substrate for the metabolism of microbial degraders. Therefore, the proper amount of amendments should be taken into account in composting to balance the motivating effect on microorganisms and the competing effect with pollutants [109]. Similarly, Hickman and Reid [110] concluded that the compost additions combined with earthworms at a ratio of 1:0.5–1:1 (soil/compost, *w*/*w*) were efficient in enhancing the removal of extractable petroleum hydrocarbons and PAHs. However, when higher ratios of compost (1:2 and 1:4) were used, PAH losses were not advanced, which may indicate that the activity of earthworms were restricted by a higher addition of compost. Wu et al. [66] showed that compost addition is an effective approach for enhancing PAH

removal from soils, especially for higher initial PAH concentrations, but increasing the ratio of compost added does not necessarily help to increase removal.

Experiments with different ratios of contaminated soil to green waste from 0.6:1 to 0.9:1 have demonstrated that in general, PAH removal is significantly enhanced in reactors increased with green waste until a maximum mixing ratio of 0.7:1 [89]. The same observation was found in Sayara et al. [78]. These results imply that low mixing ratios were not sufficient to stimulate the microbial growth; on the other hand, excessive amounts could eventually inhibit the targeted contaminants, as microbial communities prefer to use more available and easily degradable nutrients. Furthermore, co-composting of sediments (S) polluted by PAHs with urban green waste (GW) was performed using two mixing ratios (1:1 and 3:1; S/GW). In the first six months of treatment, the PAH concentrations in the 1:1 and 3:1 ratio scenarios was reduced by 57% and 26%, respectively. Despite the fact that only two mixing ratios were tested, the results again demonstrate that the low mixing ratio (3:1) was not sufficient to enhance the degradation process [94]. When different corn straw ratios (1%, 2%, 4%, or 6% *w*/*w*) were investigated for the remediation of aged PAHs in soils, removal rates were significantly (*p* < 0.05) enhanced under the 6% ratio, mainly for HMW PAHs. This indicates that the high amendment of corn straw was a potential option for the remediation of PAH-contaminated soils [111].

#### *4.5. Bioaugmentation*

When the indigenous microbial activity is not sufficient, or does not have the potential to achieve the set goals for bioremediation [112,113], it appears imperative to accelerate the process using different approaches. Among these approaches is bioaugmentation. The mechanism of this approach fundamentally depends on introducing exogenous microorganism strains that are characterized by their high capacity and diverse metabolic profiles in degrading the target contaminants [16,25,29,31,32]. However, and as concluded in many studies, the application of this approach has not always been effective in enhancing biodegradative capacity, mainly during the composting of contaminated soils [25,88,114,115]. For instance, 84% of petroleum hydrocarbon was degraded when *Candida catenulate CM1* was used as an inoculant, while a removal rate of only 48% was obtained without inoculation, indicating a positive impact of bioaugmentation [29]. On the other hand, treatments using different substrates (mixing ratio = 1:1) were performed at the laboratory and field scales, and incubated with/without fungal inoculum (*Phanerochaete velutina*). Laboratory scale treatment showed that HMW PAHs were degraded significantly in the fungal-inoculated microcosms, such that 96% of four-ring PAHs and 39% of five- and six-ring PAHs were removed in three months, whereas 55% of four-ring PAHs and only 7% of five- and six-ring PAHs were degraded in non-inoculated ones. However, the field scale achieved similar degradation rates. Importantly, the number of gram-positive, PAH-ring, hydroxylating dioxygenase genes in the field scale experiment was found to increase 1000-fold, indicating that bacterial PAH degradation played a major role [116]. Wu et al. [117] compared bioaugmentation using *Acinetobacter SZ-1* strain and biostimulation using (NH4)2SO<sup>4</sup> and KH2PO<sup>4</sup> in a petroleum-contaminated soil. It was found that the dissipation of total petroleum hydrocarbons (TPH) and the amounts of cultivable TPH, alkane, and PAH-degrading microorganisms were higher for biostimulation than for bioaugmentation. Similarly, Canet et al. [118] demonstrated that fungal inoculation, including four well-known PAH-degrading microorganisms (*P. chrysosporium* IMI 232175, *Coriolus versicolor* IMI 210866, *Pleurotus ostreatus* IMI 341687, and Wye isolate #7) in a mixture composed of non-sterile, coal-tar-contaminated soil and wheat straw, was unsuccessful in improving PAH removal. Sayara et al. [25] reported that the introduction of the white-rot fungi *T. versicolor* ATCC 42530 was not able to improve the decomposition of PAHs; on the contrary, organic substrates were capable of achieving significant degradation rates. Furthermore, inoculation with *P. chrysosporium* in a soil composting system was ineffective at enhancing the removal of benzo[a]pyrene [119].

Actually, several biotic and abiotic barriers have been documented to be behind the failure of bioaugemtation, mainly during field application of this mechanisms [51,120–122]. Biotic factors, including competition between indigenous and exogenous microorganisms for nutrients and the

biodiversity of indigenous microorganisms, could act as a barrier to the invasion of exogenous microorganisms, in addition to antagonistic interactions and predation by protozoa and bacteriophages. Abiotic factors include all the physicochemical properties of pollutants and soils, such as pH, contaminant concentration, soil type, temperature, humidity aeration, nutrient content, and redox potential.

#### **5. Bioavailability of PAHs**

In some cases, when optimal conditions for microbial degradation are provided but low or even no degradation take place, the bioavailability of the pollutant would be the most probable reason for disabling the process from proceeding forward. Actually, the bioavailability of PAHs is directly linked to the intrinsic relationship between physicochemical and microbiological factors within the composting matrix. In particular, this factor determines the fraction of the chemical compound in the soil that can be utilized or transformed by living microorganisms [68,123–126].

Bioremediation is governed by PAH sorption onto the soil matrix in such a way that gradual sorption diminishes the possibility of desorption, and thus the PAH overstates its persistency within the soil organic matrix. This would explain the biphasic behavior of contaminants during bioremediation processes, which are associated with high removal rates in the initial phase, which is primarily limited by microbial degradation kinetics; in the second phase, though, the removal rate is low and generally limited by slow desorption. PAHs with low bioavailability are characterized with low desorption mainly in the second phase of bioremediation [127,128].

#### *5.1. Factors A*ff*ecting PAH Bioavailability*

PAHs are characterized by their high hydrophobicity, consequently increasing their affinity for being adsorbed into soil organic matter and ultimately being less available for biological uptake. Different studies [70–72,129,130] have highlighted that the following factors have an essential role in determining the bioavailability of PAHs. First is contamination time (ageing): the irreversible sorption of PAH is exponentially proportional with contact time, thus decreasing the bioavailability of pollutants to microorganisms and therefore the rate and extent of biodegradation. For instance, the removal efficiency of anthracene from freshly- and age-spiked agricultural soil was investigated. The results revealed that 72% of anthracene was removed in freshly-spiked soil, while only 34% was degraded in aged soil [131]. Nevertheless, it is worth mentioning that recently contaminated soil would exhibit toxicity to or even inhibit indigenous microorganisms until they adapted to the new environment [132,133]. The second factor for determining the bioavailability of PAHs is their physicochemical properties: PAHs' water solubility is considered a crucial factor regarding their bioavailability. It is inversely proportional with PAHs' molecular weight, which in turn reduces their accessibility to microorganisms (Section 2.1). The last factor is the physicochemical properties of the soil: organic matter, particle size, and shape have a major influence on PAHs' bioavailability. Mineral surfaces (i.e., clays) and organic matter of the soil matrix are characterized by their high affinity to adsorb PAHs.

The addition of organic co-substrates to the composting mixture is believed to enhance the bioavailability of PAHs, which consequently increases the biodegradation rate [25,50,71]. Kobayashi et al. [106] demonstrated that water-extractable organic matter (WEOM) from cow manure compost was observed to increase the apparent solubility of phenanthrene, pyrene, and benzo[a]pyrene to 8.4, 34, and 89 times higher than their measured water values, respectively, thus promoting their solubility and biodegradation. Additionally, in a diesel-spiked soil, Wu et al. [66] showed that compost addition initially decreased PAH removal by up to 89% because of the decreased bioavailability resulting from strong sorption. However, as time increased, compost amendment enhanced PAH removal by more than two-fold compared with unamended soil, to which 30% was contributed by desorption and 70% by degradation. In coal tar- and coal ash-contaminated soils, compost addition

was beneficial overall for enhancing PAH removal up to 94%, and 40% of the total loss was due to the enhanced desorption [66].

The stability of the used co-substrates plays a major role in stimulating bioavailability and biodegradation of PAHS (as discussed in Section 4.3). This type of substrates was found to have more humic matter [71]. In this context, humic matter was able to increase the microbial activities much more than those developed in humin (aged organic matter), demonstrating that humin is able to sequester organic contaminants in a stronger way [70]. An important observation is that the bioavailability of the more readily degradable or LMW PAHs was decreased due to competitive inhibition of the enzymes, which is associated with biodegradation when the enzymes present in a multiple-PAH mixture. However, the bioavailability of those usually more recalcitrant PAHs (HMW PAHs) was increased by producing inducible enzymes for catabolism [134].

#### *5.2. Surfactant*

As mentioned in the previous sections, some of PAHs are characterized by their high hydrophobicity as well as low solubility, as they have the ability to be strongly adhere to soil particles and be slowly released into the water phase [135]. Among the different alternatives to overcome the problems of low bioavailability during bioremediation of PAHs is the application of surfactants. The functionality of these additives basically depends on reducing interfacial surface tension, and thus increasing their solubility [30,33,34,136,137]. The efficiency of these surfactants is influenced by many factors, including surfactant type and concentration, PAH hydrophobicity, temperature, pH, salinity, dissolved organic matter, and microbial community. An imperative and crucial element for effective PAH remediation is the selection of the optimum ratio of mixed surfactants to avoid the inhibition of microbial activities [33,34,137]. Nowadays, various groups of surfactants are available, and each one is being used under certain conditions to be compatible with the available environment. Furthermore, biosurfactants that are produced by microorganisms are receiving more favor, as they considered more environmentally friendly [33,138]. Interestingly, Both Tween-80 and rhamnolipid were found to improve the bioremediation fluoranthene [139]; however, it should be considered that the application of surfactants may not always lead to enhanced PAH biodegradation or removal. In fact, if the surfactant is preferentially used as an easier carbon substrate than PAHs for soil microorganisms, it may actually inhibit PAH biodegradation. Selection of surfactant types is therefore crucial for the effectiveness of surfactant-aided bioremediation of PAH-contaminated soils [140].

#### **6. PAH Biodegradation Pathway**

As illustrated in the literature, a wide spectrum of microorganisms has been classified, and these microorganisms are known for their high potential in degrading PAHs. These microorganisms include, but not limited to, bacteria, fungus, actinomycetes, protozoa, and algae [69,87]. Actually, the biodegradation of PAHs has the possibility of being undertaken either under aerobic or anaerobic conditions. However, aerobic conditions are more preferable, due to their documented efficiency [141]. As a result, composting as an aerobic technique has received more attention for treatment for such types of pollution [25,78,79,123,132,141].

Fortunately, it has been documented that a wide variety of bacterial cultures have the potential to biodegrade LMW PAHs directly, using them as the sole carbon and energy source [142–145]. Otherwise, PAHs (like HMW PAHs) have to proceed through the accumulation of these compounds in the body of microorganisms, and then be decomposed through sequential steps and multiple routes (Figure 2) into a bioavailable form that could be metabolized by microorganisms [123,141,146,147]. Hydroxylation of the aromatic ring via a di- or monooxygenase enzymes or dehydrogenase is the first step in the degradation process, with the formation of a *cis*-dihydrodiol, which gets rearomatized to a diol intermediate by the action of a dehydrogenase. These diol intermediates may then be cleaved by intradiol or extradiol ring-cleaving dioxygenases through either an ortho-cleavage or meta-cleavage pathway, leading to intermediates such as catechols and protocatechol acid that are ultimately converted to tricarboxylic

acid cycle intermediates, which could be considered as the end of the biodegradation [123,144,146–149]. Bacteria can also degrade PAHs via the cytochrome P450-mediated pathway, with the production of trans-dihydrodiols.

**Figure 2.** Bacterial and fungal biodegradation pathways of PAHs.

It should be noted that HMW PAH degradation pathways still need more investigation, as few bacterial isolates were found to be capable of degrading them. Also, their biodegradation in some cases is complicated and passes through different routes, or even proceeds via co-metabolism, like that of benzo[a]pyrene [150–152].

Fungal enzymatic activity also has a key role in the biodegradation of PAHs. Lignolytic and non-lignolytic fungi have the capability of oxidizing PAHs utilizing cytochrome P-450 monooxygenase and a lignin-degrading enzyme system for oxidizing aromatic rings. Usually, an oxygen atom is incorporated into the aromatic nucleus, whereas the remaining atom is reduced to water to yield cis-transdihydrodiols. The formed arene oxide, though non-enzymatic, can undergo some rearrangement to form a phenol, which can be further conjugated with glucose, xylose, gluconeric acid, and sulfate. On the other hand, ligninolytic fungi, which are usually known as white rot fungi, have been characterized by their capability to degrade PAHs through ligninolytic and non-ligninolytic culture conditions. Ligninolytic enzymes oxidizes the PAH ring by producing hydroxyl free radicals by the donation of one electron; consequently, PAH–quinones and acids are formed instead of dihydrodiols. [123,153–155]. Extracellular enzymes of white rot fungi, which include laccase, LiP, and MnP, have a key role in the degradation of PAHs [153,156–159].

#### **7. Conclusions**

The main conclusion of this review is that the use of compost and composting in several strategies significantly improves the removal of PAHs in contaminated soils. However, this strategy should be well studied and tested. For instance, future studies are required on compost stability, as it is an important parameter for considering the removal of PAHs. Composting also needs to be optimized to improve PAH removal. This could include, for instance, the use of some additives like surfactants, which can be of help for the desorption and further removal of PAHs. Furthermore, more investigations are still needed regarding the biodegradation of PAHs combined with other hydrocarbons in mixtures, biodegradation HMW PAHs, and the microbial interactions within PAH-degrading consortia. In summary, composting and compost opens a wide number of strategies to improve the bioremediation of PAH-contaminated soils. However, it is important to define this strategy and to test its efficiency before full-scale application.

**Author Contributions:** Writing—original draft preparation, T.S.; writing—review and editing, A.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** Tahseen Sayara would like to thank Palestine Technical University for its administrative support.

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

#### **Abbreviations**


#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Chemical and Spectroscopic Investigation of Di**ff**erent Soil Fractions as A**ff**ected by Soil Management**

#### **Francesco De Mastro, Claudio Cocozza, Gennaro Brunetti \* and Andreina Traversa**

Dipartimento di Scienze del Suolo, della Pianta e degli Alimenti, University of Bari, Via Amendola 165/A, 70126 Bari, Italy; francesco.demastro@uniba.it (F.D.M.); claudio.cocozza@uniba.it (C.C.); andreina.traversa@uniba.it (A.T.)

**\*** Correspondence: gennaro.brunetti@uniba.it

Received: 14 March 2020; Accepted: 5 April 2020; Published: 9 April 2020

**Abstract:** The interaction of organic matter with the finest soil fractions (<20 µm) represents a good way for its stabilization. This study investigates the effects of conventional (CT), minimum (MT), and no (NT) tillage, fertilization, and non-fertilization, and soil depth (0–30, 30–60, and 60–90 cm) on the amount of organic carbon (OC) in four soil fractions. Diffuse reflectance infrared Fourier transform spectroscopy (DRIFT) was performed to obtain information about the OC quality and the mineralogical composition of these fractions. The CT shows the highest amount of the finest fraction while the fertilization enhances the microbial community with the increase of soil micro-aggregates (250–53 µm). The coarse fraction (>250 µm) is highest in the upper soil layer, while the finest fraction is in the deepest one. The greatest OC content is observed in the topsoil layer and in the finest soil fraction. DRIFT analysis suggests that organic components are more present in the finest fraction, calcite is mainly localized in the coarse fraction, quartz is in micro-aggregates and 53–20 µm fraction, and clay minerals are in the finest fraction.

**Keywords:** tillage; fertilization; soil depth; organic carbon; clay minerals; diffuse reflectance; infrared Fourier transform spectroscopy

#### **1. Introduction**

Previous studies have used the particle size fractionation for obtaining information about the influence of land use and depth on the distribution of soil organic carbon (SOC) [1,2]. The various soil fractions can differently immobilize organic carbon (OC) through the formation of organo-mineral complexes [3,4]. In particular, quartz particles exhibit only weak bonding affinities to organic matter (OM), while clay size particles (i.e., sesquioxides and phyllosilicates) have a large surface area and numerous sorption sites [3,4]. The physical protection of OM through its occlusion within clay minerals limits its microbial decomposition, which reduces the C mineralization [5]. Therefore, the sand related OM represents the active pool of soil organic matter (SOM), the OM linked to the silt size fraction is the intermediate pool, and the clay related OM represents the passive and the older SOM pool [4]. The microaggregates, composed mainly of clay minerals, represents the most efficient way to stabilize the SOM [6,7] by forming bridges between the exchangeable cations of layer silicates and functional groups of organic compounds [8]. The formation of macroaggregates is favored by the decomposition of fresh plant residues and fungal hyphae [9]. The SOM in the macroaggregates is available for microbial utilization while the protected microaggregates form a long-term reserve of mineral-associated C that is not "humified" and can be attacked by microorganisms once exposed [9].

It is known that chemical fertilizers can modify soil physical, chemical, and biological properties with clear consequences on soil aggregates [10]. A different development of the root system, stimulated by fertilization, influences the production and release of root exudates, directly involved in the formation of aggregates [11]. In addition, a greater growth of the root system following the fertilization increases the SOM [11]. Several previous studies have demonstrated that tillage practices influence the content and the dynamic of SOM [7,12–14]. Soil tillage increases the turnover of macroaggregates by inhibiting the formation of microaggregates within macroaggregates and, thus, reducing the sites where the OM is stabilized [5].

In order to identify the minerals and their changes among the different soil fractions, the diffuse reflectance infrared Fourier transform (DRIFT) spectroscopy has been utilized. This technique can be considered rapid, inexpensive, and precise, and can be applied to estimate the water-bearing minerals, such as clay minerals together to other sheet silicates as muscovite, illite, smectite, kaolinite, and chlorite [15], and the presence of organic matter [16]. DRIFT spectroscopy allows us to analyze the matrices without pressing them by avoiding the error due to scattering [17], and to have a band intensity four times greater than that of IR spectroscopy due to the non-mixing of the soil sample with KBr [18].

The objective of this study was to investigate the effects of different soil managements on the quantity of SOC associated with several soil-size fractions. In addition, with DRIFT analysis, we tried to better understand the OM interaction with the mineral parts of the different soil fractions.

#### **2. Materials and Methods**

#### *2.1. Study Area and Experimental Design*

The trial was conducted in the experimental station of the University of Bari (Italy) located at Policoro (40◦10′20" N; 16◦39′04" E; altitude: 15 m above sea level). The soil texture was classified as silt loam (sand 8%, silt 68%, clay 24%), according to the USDA [19]. Since 2005, a two-year rotation of durum wheat with faba bean in a split-block design with three field replications has been introduced. Treatments were as follows: i) no tillage and no fertilization (NT), ii) NT and crop fertilization (30 kg P2O<sup>5</sup> ha−<sup>1</sup> : NTF), iii) minimum tillage (20 cm deep subsoiling in late August and 15 cm deep disk harrowing in November) and no fertilization (MT), iv) MT and crop fertilization (MTF), v) conventional tillage (35 cm deep moldboard plowing in late August and 15 cm deep disk harrowing in November) and no fertilization (CT), and vi) CT and crop fertilization (CTF). More details about treatments and soil properties are reported elsewhere [20,21].

After more than a decade of trial, each faba bean plot (30 × 30 m) was sampled in July 2017 at three different depths (0–30, 30–60, and 60–90 cm) using an auger, after the harvest, and the removal of aboveground crop residues. Due to the soil homogeneity, nine sub-samples have been collected from each plot using a grid sampling scheme.

#### *2.2. Particle-Size Physical Fractionation and Determination of Organic Carbon in Fractions*

Fraction-size separation was obtained by ultrasonic dispersion, according to Amelung and Zech [22], and wet sieving, according to Bornemann et al. [23] (Figure 1).

About 150 mL of milli-Q® ultrapure water were added to 30 g of air-dried soil and the suspension was gently sonicated by placing the probe tip 15 mm below the water surface and using a probe-type sonicator Sigma Aldrich, 500-Watt model (60 J mL−<sup>1</sup> ). This weak sonication was used for preserving micro-aggregates from disruption [24]. The first fraction (macro-aggregates fraction, A: 2000–250 µm) was separated from the suspension by wet sieving (250 µm), and the filtered remnant was sonicated a second time at 440 J mL−<sup>1</sup> and separated by wet sieving using sieves with different meshes (53 µm and 20 µm). After this step, the obtained fractions were: fraction B (microaggregates fraction, 250–53 µm), fraction C (coarse silt-sized fraction 53–20 µm), and fraction D (free fine silt plus clay fraction, <20 µm). All fractions were dried at 35 ◦C before elemental analysis. The water-soluble organic fraction was

isolated and discharged from each D fraction, according to Zsolnay [25], to avoid any interference. Briefly, an aliquot of each air-dried D fraction was suspended in water (1:10, w/v), and mechanically shaken for 15 min. The suspensions were then centrifuged at 6000 rpm for 15 min and the supernatant was removed.

**Figure 1.** Fractionation scheme of soil.

− The OC content of all soil fractions was determined in triplicate using a Flash 2000 CHNS-O Elemental Analyser (Thermo Scientific) calibrated by an organic analytical standard consisting of a low organic content soil with 1.55% (w/w) of carbon. About 6–7 mg of each soil fraction were dried at 40 ◦C and pre-treated with hydrochloric acid (HCl 1%) to dissolve carbonates. The OC stocks were calculated by multiplying the C concentrations and the corresponding particle-size weights.

μ

#### *2.3. Spectroscopic Analysis of Particle-Size Fractions*

μ

<sup>−</sup> μ μ μ μ Diffuse Reflectance Fourier Transform (DRIFT) spectra were recorded for each fraction in triplicate and in transmittance mode using a Thermo Nicolet Nexus FT-IR spectrophotometer, which was equipped with a Nicolet Omnic 6.0 software. Before DRIFT analysis, air-dried samples were thoroughly

μ

mixed to obtain a representative sample and then finely ground in a mill. About 200 mg of the mixture was filled in a cup and the surface was smoothed with a plastic slide. Spectra were recorded in the range of 4000 to 400 cm−<sup>1</sup> , with 4 cm−<sup>1</sup> resolution and 16 scans min−<sup>1</sup> for each acquisition.

#### *2.4. Statistical Analysis*

All analyses performed on soil fractions were conducted in triplicate. The analysis of variance (four-way ANOVA) and the Tukey's test (R software, version 3.2.3) were used to measure the effect of fertilization, tillage, and depth on the OC content for each soil fraction.

#### **3. Results and Discussion**

#### *3.1. E*ff*ects of Treatments on the Amount of Each Soil Fraction and on their Organic Carbon Content*

The physical fractionation recovered 98% of the mass and 99% of the OC from all samples. Such percentages were comparable to those previously reported by other authors [26–28], which indicates that the loss of material was very low and confirms the efficiency of the fractionation method adopted.

Table 1 shows the amounts of soil dry matter obtained from each fraction (g kg−<sup>1</sup> soil), while Table 2 reports the corresponding statistical analyses, as influenced by treatments (soil depth, tillage, and fertilization).



CT: Conventional tillage. CTF: Conventional tillage fertilized. MT: Minimum tillage. MTF: Minimum tillage fertilized. NT: No tillage. NTF: No tillage fertilized.

On average, 70% to 75% of the soil fractions consisted of fine silt and clay (<20 µm). The A, B, and C fractions decreased with depth (Table 1), while the D fraction showed an inverse trend with the only exception being the deepest layer of NTF treatment. The quantity of the smallest fraction (<20 µm) increased with depth and ranged from 500 to 867 mg kg−<sup>1</sup> .

The fertilized plots resulted in the highest amount of B fraction and the lowest amount of D fraction with respect to the unfertilized ones. Since microaggregates are the result of microbial decomposition of SOM from the macroaggregates [9,29], the highest amount of B fraction in fertilized soils could derive from the higher microbial activity promoted by the same fertilization. For example, Liao et al. [30] found the higher fungal abundance in microaggregates (250–53 µm) regardless of the type of fertilization.


**Table 2.** Analysis of variance and mean values of the amount of each soil fraction subdivided by soil depth, tillage, and fertilization treatment. The standard deviation is in parentheses.

CT: Conventional tillage. MT: Minimum tillage. NT: No tillage. The values in each column followed by a different letter are significantly different according to Tukey's test. n.s.: not significant. \*\*\* significant at the *P* ≤ 0.001.

With regard to the tillage, the highest amount of D fraction and the lowest amount of B fractions were obtained from CT soils, which is likely due to the major physical disturbance and microbiological activity induced by the conventional tillage that increased macro-aggregates and micro-aggregates turnover [5].

Table 3 shows the analysis of variance and mean values of OC content of soil fractions, as affected by soil depth, tillage, fertilization, and size fraction. The interactions among these parameters were not significant (data not shown) except the one between the OC content of each soil fraction and soil depth (*P* ≤ 0.001) since, as expected, the OC content of all fractions decreased with depth.

The significant decrease of the OC content of all fractions from the upper layer (on average, 125.9 mg OC kg−<sup>1</sup> fraction) to the deepest one (54.5 mg OC kg−<sup>1</sup> fraction) resembles the common stratification of SOC along the profiles. The highest value of OC was found in the D fraction due to the entrapment of organic components in the finest fractions of soil, as reported by Gregorich et al. [31]. This highlighted the role of clay particles in the OM stabilization due to their high specific surface area and charge. In fact, clay minerals are considered the most active constituents in the formation of organo-mineral complexes [32] and are responsible for long-term preservation of soil OM, even over millennia [33].

No kind of tillage significantly affected the OC content of each fraction likely due to the balance in the soil achieved because of the long-term experiment, as reported by Rita et al. [34] who found no significant difference in OC fraction content among several 30-year-old land-use trials. In addition, the main OM input, which is the above-ground crop residues, has been removed from the field at the end of the crop cycles in all treatments. With regard to the latter topic, it has been demonstrated that the aboveground crop residues are important for building up soil fertility not only as input of OM, but also because they cover the soil during the hot weather, which conserves SOM [35].


**Table 3.** Analysis of variance and mean values of the OC in soil fractions, subdivided by soil depth, tillage, fertilization treatment, and size of fractions. The standard deviation is in parentheses.

CT: Conventional tillage. MT: Minimum tillage. NT: No tillage. A (>250 µm). B (250–53 µm). C (53–20 µm). D (<20 µm). The values in each column followed by a different letter are significantly different according to Tukey's test. n.s.: not significant. \*\*\* Significant at the *P* ≤ 0.001.

Lastly, the fractional OC content did not differ between fertilized and not fertilized plots since the adopted fertilization was only inorganic and rather low.

#### *3.2. E*ff*ects of Treatments on the Spectroscopic Properties of Each Soil Fraction*

Figure 2 shows the DRIFT spectra of different soil fractions under various tillage.

The peak at about 3623 cm−<sup>1</sup> can be ascribed to the O-H vibration in the octahedral layers of 2:1 and/or 1:1 silicates [36]. The same peak can be highlighted by removing the SOM through an appropriate procedure [37]. This peak was observed in B, C, and D fractions regardless of the type of tillage, with a slight increase of the relative intensity as the fraction size decreased. This was in line with other papers [38–41] showing that phyllosilicates (kaolinite, chlorite, smectite, illite) are the main components of the clay fraction of soils. The peak at 2927 cm−<sup>1</sup> , ascribed to the stretching of the aliphatic C-H group, was evident only in the A fraction, and could be due to the signal of organic matter consisting of labile plant residues [42]. Additionally, in this case, this peak can be highlighted by removing minerals through HF treatment [43]. The peak at about 2517 cm−<sup>1</sup> can be attributed to the CO<sup>3</sup> stretching and calcite bending, as indicated by peaks at about 1450, 867, and 698 cm−<sup>1</sup> [15,44]. These peaks, especially those at 1450 and 867 cm−<sup>1</sup> , were more pronounced in the spectrum of the A fraction. High percentage of calcite in the sand-sized fraction of a Mediterranean soil was also found in a previous work [3] and was ascribed to its lithogenic origin. The peaks at about 1991, 1868, 1793, and 698 cm−<sup>1</sup> were related to Si-O bending of quartz minerals [15,43,44] as well as the peaks at about 791 cm−<sup>1</sup> related to the Si-O stretching of quartz minerals [15,43]. Overall, these peaks were slightly more pronounced in all fractions B and C possibly because of the physical breakup of sand size quartz into silt dimension. In contrast, the reduced intensity of the same peaks in D fractions suggested a limited physical alteration of quartz minerals in the finest particles and, therefore, an intermediate stage of soil evolution [45]. The peak at about 1630 cm−<sup>1</sup> was ascribed to aromatic C=C skeletal vibrations, C=O stretching of quinone and amide groups, C=O of H-bonded conjugated ketones, and it was typical of the organic components [46]. As expected, this peak was absent in the A fraction, appeared in the B fraction, and became more pronounced from fraction C to D. The presence of the peak linked to aromatic structures in the fraction D could be due to the presence of organic matter involved in the formation of organo-mineral associations and in the coating of the mineral surface by sorption or precipitation processes [47–49]. The broad band at about 1030 cm−<sup>1</sup> can be ascribed to the stretching of the carbohydrate and polysaccharides-like substances [50]. It was found mainly in the silt-clay and free fine silt plus clay fractions (fraction C and D). Many previous studies have reported high proportions of these compounds in the mineral-associated organic matter fraction [5,51–54]. Polysaccharides of microbial origin mainly bind clay particles by promoting the formation of microaggregates of <50 µm [55]. Glicoproteins of fungal origin, such as glomalin, contain about 85% of sugars and are decomposed very slowly in soil [56]. In contrast, the lower relative intensity of the previously mentioned peak in fraction A suggested the presence of polysaccharides of plant origin responsible for the formation of easily degradable macro-aggregates [56]. The bands at about 529 and 478 cm−<sup>1</sup> can be ascribed to Si-O-Al and Si-O-Si vibrations, respectively, and are distinctive of phyllosilicates [41]. The relative intensity of the second band decreased in proportion to the size of the fractions. Ndzana et al. [41] reported similar results suggesting that the crystalline structure of phyllosilicates weakened in the finest soil fraction.

DRIFT spectra recorded from NT, MT, and CT samples were similar. Therefore, Figure 3 shows only the DRIFT spectra of each soil-size fraction isolated from the NT treatment along the soil profile. The only slight difference evident in all the soil-size fractions was the peak among 1493 and 1450 cm−<sup>1</sup> , which increased its relative intensity with depth. This suggests a slight increase of calcite along the soil profile due to dissolution/precipitation phenomena typical of aridic climates and high soil pH [45]. In addition, a slight reduction of the relative intensity of the peak at 1628 cm−<sup>1</sup> was observed only in the D fraction, according to the reduction of organic matter content with depth.

Figure 4 reports the DRIFT spectra of each soil fraction isolated from the 0–30 layer of the NT treatment fertilized and unfertilized. The spectra of fractions A, B, and C were very similar between the two levels of fertilization. The fraction D of NTF treatment showed a greater relative intensity of the peak at 1027 cm−<sup>1</sup> compared to NT treatment, possibly due to the greater quantity of polysaccharides coming from microorganisms whose activity is certainly favored by fertilization, as reported by De Mastro et al. a [20].

Figure 5 shows the DRIFT spectra of each soil fraction isolated from the 0–30 layer under NT, MT, and CT practices. Even in this case, the spectra of fractions A, B, and C were more similar regardless of the kind of tillage. The peak at 1027 cm−<sup>1</sup> showed a greater relative intensity in fraction D of MT and CT treatments compared to NT since the higher aeration of the formers enhanced the microbial activity.

**Figure 2.** Diffuse reflectance infrared Fourier transform spectra of all soil fractions isolated from different tillage treatments at 0–30 cm of depth.

**Figure 3.**Diffuse reflectance infrared Fourier transform spectra of all soil fractions (A, B, C, and D) isolated from no tillage (NT) treatment at different soil depths.

**Figure 4.** Diffuse reflectance infrared Fourier transform spectra of the different soil fractions (A, B, C, and D) isolated from no tillage (NT) treatment (0–30 cm) fertilized and not fertilized.

**Figure 5.** Diffuse reflectance infrared Fourier transform spectra of all soil fractions (A, B, C, and D) isolated from different tillage treatments at 0–30 cm of depth.

#### **4. Conclusions**

The soil management influenced the quantity of soil fractions since the CT enhanced the finest ones, whereas the fertilization increased the B fractions, which possibly fuels the development of a microbial community that fosters microaggregate formation. In general, soil depth influenced the amount of each fraction, with higher amounts of fraction A in the upper soil layer and higher amounts of the finest one is more evident in the deepest. The OC content was primarily influenced by fraction size and soil depth since higher OC content was found in the topsoil layer (0–30 cm) and in the finest soil fraction (fraction D), as confirmed by the DRIFT analysis. The different tillage may increase the mass of the soil fractions but not their OC content. However, MT and CT affected positively the quality of the OM stimulating a relative major production of polysaccharides of microbial origin that possibly stabilize the finest size fractions due to their major recalcitrance. The same result is obtained with fertilization. The DRIFT analysis can provide information about the quality of the main minerals present in the different soil-size fractions, since fraction A appeared mostly rich in calcite, fractions B and C appeared mostly rich in quartz, and the finest fraction showed the highest content of phyllosilicates. The different quality of minerals among the soil fractions suggested an early-intermediate stage of weathering and did not change with tillage and fertilization, except for calcite whose relative intensity increased with depth.

**Author Contributions:** Conceptualization, G.B., F.D.M., A.T., and C.C. Methodology, G.B., F.D.M., A.T., and C.C. Software, G.B., F.D.M., A.T., and C.C. Validation, G.B., F.D.M., A.T., and C.C. Formalanalysis, F.D.M. Investigation, F.D.M. Resources, G.B. and C.C. Data curation, G.B. and A.T. Writing—original draft preparation, A.T. Writing—review and editing, G.B., F.D.M., A.T., and C.C. Visualization, G.B. and C.C. Supervision, G.B. Project administration, G.B. Funding acquisition, G.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

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

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
