Going Beyond Soil Conservation with the Use of Cover Crops in Mediterranean Sloping Olive Orchards

: Among the agricultural practices promoted by the Common Agricultural Policy to increase soil functions, the use of cover crops is a recommended tool to improve the sustainability of Mediterranean woody crops such as olive orchards. However, there is a broad range of cover crop typologies in relation to its implementation, control and species composition. In that sense, the inﬂuence of different plant species on soil quality indicators in olive orchards remains unknown yet. This study describes the effects of four treatments based on the implementation of different ground covers (CC-GRA: sown cover crop with gramineous, CC-MIX: sown cover crop with a mixture of species and CC-NAT: cover crop with spontaneous vegetation) and conventional tillage (TILL) on soil erosion, soil physicochemical and biological properties after 8 years of cover crop establishment. Our results demonstrated that the presence of a temporary cover crop (CC), compared to a soil under tillage (TILL), can reduce soil losses and maintain good soil physicochemical properties and modify greatly the structure and diversity of soil bacterial communities and its functioning. The presence of a homogeneous CC of gramineous ( Lolium rigidum or Lolilum multiﬂorum ) (CC-GR) for 8 years increased the functional properties of the soil as compared to TILL; although the most relevant change was a modiﬁcation on the bacterial community composition that was clearly different from the rest of treatments. On the other hand, the use of a mixture of plant species (CC-MIX) as a CC for only two years although did not modify greatly the structure and diversity of soil bacterial communities compared to the TILL soil, induced signiﬁcant changes on the functional properties of the soil and reverted those properties to a level similar to that of an undisturbed soil that had maintained a natural cover of spontaneous vegetation for decades (CC-NAT).


Introduction
Currently, soil degradation is caused by several processes, of which erosion has been identified as the most serious threat to soil health maintenance [1]. Diverse ecosystem services and soil functions associated with this natural resource are in danger as current soil loss rates are orders of magnitude higher than natural soil formation [2]. These rates have a large variability under the same land use, since erosion processes are highly affected by a combination of specific environmental factors such as climate, soil type and topography and management operations [3].
Soil managements based on tillage or the use of herbicides to control adventitious vegetation has been proved unsustainable in terms of soil conservation, since it can reach soil losses systematically above 20 t·ha −1 ·year −1 on Mediterranean sloping areas in the short [4] and long term [5]. Cover crops have been proven as an alternative to conventional tillage in woody crops, such as olive, for soil management. Thus, currently, there is a tendency in sloping (slope > 10%) olive production systems, mainly induced by the Common Agricultural Policy obligations, to reduce bare soil management through the implementation of permanent or temporary cover crops [6,7]. The implementation of cover crops has also been promoted as a way to enhance associated ecosystem services [8,9] and are highly recommended as a soil and water conservation strategy for their agronomical and environmental benefits. Several studies have demonstrated multiple benefits, in the long term, of using temporary, from fall to early spring, cover crops in olive orchards at Mediterranean sloping areas, including an increase in carbon sequestration, soil structure improvement (soil macroporosity, soil aggregates and infiltration rates), increase in topsoil organic matter, nitrogen and macronutrients content (e.g., available P and K), reduction of compaction and water erosion control [4,[10][11][12][13][14][15][16]. In addition, a favorable correlation between ground cover maintenance under proper control to avoid the risk of competition with the tree has been shown, with an increase on tree growth and yield and physiological performance, being the treatments that enhance nutrient absorption by olives trees, the most profitable ones [17][18][19].
In the last decade, the use of biological indicators has gained importance for assessing soil health (e.g., micro-meso fauna, enzymatic activities, microbial biomass, etc.). As a contrast to the relatively high number of studies describing the benefits of using cover crops on physical and chemical soil quality indicators in Mediterranean olive orchards, those assessing specifically the effects of ground cover management (also considering different plant species composition) on soil biological properties, and more specifically on the microbial community structure and its activity are less numerous, e.g., [20][21][22][23][24][25][26][27][28][29]. In all these studies, a positive effect of the use of cover crops as compared to tillage is described for some biological indicators (soil respiration and catabolic activity, assessment of specific groups of organisms, etc.), observing a general improvement of soil biological functionality. However, there is still little information on the effect of temporary and/or permanently covered managed systems in olive orchards on total soil microbial community composition and its relation to its biological activity, particularly among different types of cover crops. Moreover, the establishment of links among those biological indicators with other physicochemical ones is generally lacking.
One key aspect when implementing suitable soil management based on cover crops in woody orchards is the selection of plant species (both in mono and multispecies cover crops and either sown or spontaneous) due to several remaining uncertainties regarding the potential benefits of each species, and the different agricultural practices needed to control the vegetation [30,31]. Only a few studies have assessed that issue in detail but paying most of the attention to the physicochemical soil properties, erosion control or olive yields [10,[32][33][34].
Therefore, starting from the hypothesis that cover crops promotes and preserve soil health, the facts described above, but also that the composition of the cover crops is a major driver of soil functions. As such, the main objectives of our study were to assess whether: (i) the presence of a cover crop and its composition, regarding plant species, can significantly influence soil physicochemical and biological properties in an olive orchard under semiarid conditions; (ii) bacterial community composition in soil is modified due to the cover crop presence or its composition; (iii) those changes can be associated to some physicochemical or specific microbial activities induced in the soil due to differences on erosion processes.

Study Site
The field trial was carried out at a commercial olive orchard established in 1985 and located in Benacazón (SW Spain; 37 • 20 36" N, 6 • 13 45" W) with an average elevation of 92 m.a.s.l. ( Figure S1). Olive trees from the variety Gordal were spaced at 8 m × 6 m and were drip irrigated during the dry season since their plantation. The climate is the Mediterranean with an average annual precipitation of 534 mm, concentrated mostly in late fall and winter, and an average annual air temperature of 18.6 • C. The soil belongs to the Petrocalcic Palexeralf [35] and is well drained, with an average organic matter content of 1.3%, 28% of CaCO 3 and sandy loam texture class [4].

Soil Managements and Agronomic Practices
The study site belongs to a medium-term erosion trial in the olive orchard, originally established in 2003, conformed by six runoff plots, and aimed to compare the effect of different ground cover managements on soil quality and erosion at hillslope scale. The experimental plots are north-west oriented with an average steepness of 11% [4,36].
One of the ground managements was conventional tillage (TILL) consisting of regular chisel plow passes (2-3 times per year at 10-15 cm depth) depending on adventitious vegetation growth. At the sampling time, this plot had been for 8 years under this management.
The second management (CC-GRA) maintained a temporary cover crop (CC) of gramineous during the same period. Seeds of Lolium rigidum or Lolium multiflorum, depending on availability, were sown during early fall at 100 kg of seeds per ha along the inter-tree rows. This cover crop was chemically controlled using contact herbicide (glyphosate 36% brand name Touchdown Premium ® ) during April, with the exact date depending on the rainfall amount during late winter and early spring of the particular year. A central strip, approximately 1 m wide, was left unsprayed to allow self-seeding of the cover crop for the next year.
The third management included a soil that had been under different soil managements during the 8-year period: TILL for 2 years, CC-GRA for 4 years and a temporary cover crop for 2 years sowed with a mix of plant species (CC-MIX) including Coriandrum sativum L., Borago officinalis L., Vicia sativa L., Nigella papillosa G., Calendula arvensis L., Echium creticum L., Melilotus indicus L., Diplotaxis virgata Cav., Silene vulgaris Moench and Salvia verbenaca L. The purpose of this seed mix was to provide an expanded flowering period to enhance arthropod biodiversity [11]. This cover crop was mowed, to prevent competition for soil water with the olive trees, in mid-May.
Additionally, a fourth management (CC-NAT) parallel and nearby to the runoff plots, and within the same olive orchard, was delimited. This area has been under the continuous presence of spontaneous vegetation for more than 20 years and was sampled to serve as a reference site of a low disturbed soil within the orchard. Vegetation was controlled through regular mowing, 1-3 times a year depending on rainfall.
For the cover crop treatments, except CC-NAT, the inter-rows were fertilized at the time of seeding with nitrogen (50 kg·ha −1 ) during the fall by direct application, and additional fertilization spaced 1-3 years depending on the cover crop status. Adventitious vegetation was chemically controlled (Fluroxipir and Flazasulfuron) along the tree rows, periodically during late winter-early spring to avoid competition for water with the olive trees. Olive trees were ferti-irrigated with ammonium sulfate in February (1.6 kg·tree −1 ) and from April to October, during the irrigation season. The fertilization rates were the ones required for replenishment, established by periodic foliar analysis. The average amount of water applied during the irrigation season was 240 mm, although it slightly differed annually depending on the season's rainfall.

Erosion Measures
Runoff and sediment generated on each plot were collected after every rainfall event by three serial fiberglass tanks with flow splitters since 2003 as described in [4]. At the time of sampling, erosion measurements had been recorded for 8 consecutive years for TILL and CC-GRA treatments, and for 2 years for the CC-MIX treatment. On the first working day after a single large rainfall event, or after several small events, the collection tanks were sampled, and the runoff volume and wet sediment weight were recorded. More details on the collection system appears in [4].

Soil Sampling
Soil samples were collected during summer of the 8th year at the erosion plots, each one with a different management (TILL, CC-GRA and CC-MIX). A nearby row that had kept the spontaneous ground vegetation (CC-NAT) was sampled as well ( Figure S1). Soil samples were taken at the top 10 cm at six random points along the inter tree row slope at three slope points (upper, U, middle, M, and lower, L). Each of the six samples was mixed to produce a composite sample in each slope area sampled (U, M and L), placed in cold storage and taken back to the laboratory. Soil samples were air dried, crumbled and passed through a 2-mm sieve for the determination of physicochemical and biological properties. Soil samples were stored at 4 • C until analysis.

Bacterial Population Densities
Soil suspensions were prepared by mixing 5 g from each sample with 50 mL of sterile saline solution (0.85% NaCl, w/v). The solution was shaken at 250 rpm for 30 min using an orbital shaker and centrifuged at 1000× g for 5 min to precipitate most soil particles. This soil supernatant was used directly to estimate bacterial densities, enzymatic activities and for DNA extraction or was 1/10 diluted in 0.85% NaCl for the community-level physiological profile analysis as described below.
Culturable bacterial populations were evaluated using the dilution plating method. Tenfold serial dilutions of the soil supernatant were diluted in sterile distilled water; then, two aliquots of 100 µL of each dilution were plated on nutrient agar (Difco, Detroit, MI, USA) supplemented with cycloheximide (100 µg/mL). Plates were incubated in darkness for 5 days at 28 • C to determine the number of culturable bacteria in each soil sample.

Metagenomic Analysis of Soil Bacterial Communities
Twenty-five microliters of the soil supernatant obtained for the bacterial population density estimation was centrifuged at 12,000× g for 30 min; then, the supernatant was removed and the soil pellet containing the finest particles and microbial cells was used to extract DNA. Total DNA was extracted from three replications for each treatment using the MoBio Ultraclean™ soil DNA isolation kit (MoBio laboratories, Inc.; Carlsbad, CA, USA) according to manufacturer instructions and the FastPrep-24 Instrument (MP Biomedicals, Inc., Illkirch Cedex, France) run at 6.0 m/s for 40 s. Extracted DNA was eluted in ultrapure sterile water (Sigma, Madrid, Spain) and quantified in triplicate using the Quant-iT DNA Assay Kit Broad Range fluorometric assay (Molecular Probes Inc., Leiden, The Netherlands). DNA concentration from soil samples ranged from 20 to 60 ng/µL and was used directly for amplification.
Amplification of the V2-V3 hypervariable region of the bacterial 16S rRNA gene was performed in triplicate per DNA sample using a two-step PCR protocol with the 16S rRNA primers 27F and 357R linked to universal M13/pUC forward and reverse primers (M13F-8F and M13R-357R) [38]. Then, 10-bp barcode sequences were inserted between the full sequencing adapter A and M13F primer sequences in a second PCR reaction using a 10× dilution of the first PCR product and the reverse primer containing the full sequencing adapter B and the M13R sequence. A total of 12 barcodes were used to differentiate each of the 4 soil managements and the three slope areas sampled within each treatment. All PCR reactions, including positive and negative controls, were run in a T100TM Thermal Cycler (Bio-rad, Madrid, Spain) using the FastStart High Fidelity Polymerase (Roche Diagnostics GmBH, Mannheim, Germany) and conditions recommended by the manufacturer for pyrosequencing analysis. The PCR products were purified using the Amicon Ultra 0.5 mL 30 k, (Merck Millipore, Merck KGaA, Darmstadt, Germany), further purified twice with AgencourtH AMPureH XP PCR purification system (Agencourt Bioscience Co., Beverly, MA, USA) and quantified using the QuantiT dsDNA BR assay kit (Invitrogen, Carlsbad, CA, USA) and a fluorometer (BioTek Instruments, Winooski, VT, USA). Subsequently, all samples from each purified library were pooled in equimolar concentrations. Pools of the 12 samples were diluted to obtain a total of 1 × 10 5 copies/µL (1.0 molecules per bead), and an emulsion PCR was performed with the Lib-L kit (454 Life Sciences). DNA positive beads were enriched, counted on the GS Junior Bead Counter, and loaded onto a picotiter plate for pyrosequencing on the 454 Life Sciences (Roche, Mannheim, Germany) Junior platform according to standard 454 platform protocols. Two independent pyrosequencing runs were performed for the same library.

Soil Microbial Functional Activities
Soil respiration measurements were made with the OxiTop Control method (WTW, Weilheim, Germany) using the procedure mentioned in [37]. Measurements were based on pressure variation in a closed system due to oxygen consumption and CO 2 released by microorganisms in the sample. Briefly, 100 g of air-dried soil were moistened to 50% of soil water holding capacity and was placed at the bottom of a hermetically sealed OxiTop crystal container with a 50 mL beaker filled with a 1 M NaOH solution. The container was incubated in a dark cabinet at 25 • C for 8 days.
Nineteen soil enzymatic activities were assayed using the API ZYM TM (BioMerieux, Madrid, Spain) kit as described in [39]. Briefly, an aliquot (65 µL) of the soil supernatant was dispensed into each of the 20 microcapsules. The API ZYM strips were then covered and incubated at 28 • C for 24 h in the dark. After incubation, the galleries were activated by adding 30 µL of each reagent (ZYM A and ZYM B; BioMerieux, Madrid, Spain). After 5 min at room temperature, semiquantitative evaluation of the activities was measured referring to a colorimetric standard table by assigning a numerical value of 0-5 (equivalent to 0-40 nmol) depending on the chromogenic substrate intensity produced by the hydrolase reaction [40]. For the purposes of this study, the results were also reported as the total enzymatic activity (TEA) measured in nanometers that could be achieved.
Intensity and diversity of bacterial metabolism were evaluated as described in [39] using a community-level physiological profile (CLPP). Soil suspensions in 0.85% NaCL (150 µL) from each sample were added to each well of a Biolog EcoPlate TM system (Biolog Inc., Hayward, CA, USA) with the three independent samples taken at the different slope positions being assayed in the same microplate. Plates were incubated at 28 • C and color formation was measured at 590 nm with a Tecan Safire fluorospectrometer (Tecan Spain, Barcelona, Spain). Readings were made at regular time intervals for 192 h. Average well color development (AWCD) of each plate was calculated as the mean of the absorbance values for all 31 response wells per reading time and kinetics of AWCD was used to determine the speed and the level of development of the bacterial communities. The metabolic diversity of the bacterial communities was assessed by comparing the CLPPs of all replicates at the end of the experiment (CLPPs-Final). Normalization for the effect of putative variable inoculum densities among plates was achieved by dividing individual well responses by the AWCD of the corresponding microtiter plate before statistical analysis [39,41].

Statistical and Bioinformatic Analyses
Different statistical analyses were performed depending on the data to be analyzed. Erosion measurements (runoff and soil losses) were analyzed by ANOVA followed by a Tukey's means multiple comparison (p = 0.05) using InfoStat (Córdoba, Argentina). Data of soil physicochemical properties, soil respiration, API ZYM activities and CLPPs-216 were subjected to cluster analysis using the Bray-Curtis distance measure. Dendrograms were constructed with the Ward Method using Bionumerics 7.6 (Applied Maths, Sint-Martens-Latem, Belgium). Physicochemical soil properties, soil respiration, richness diversity statistics of metagenomics data and richness of a total number of distinct enzymatic activities (API ZYM) or carbon source assimilated profiles (CLPP) (S API ZYM and S CLPP ), respectively, were analyzed using the rank-based Kruskal-Wallis test to determine significant differences (p = 0.05) among the different soil management. When appropriate, an all-pairwise multiple comparison test on ranked data was used to determine where those differences among treatments occurred.
For bioinformatics analysis of metagenomic data, the Fastq files were checked for quality with FASTQC [42] and filtered with Cutadapt software [43], removing poor-quality, short-length reads and sequences adapters. Trimmed reads were analyzed within Quantitative Insights into Microbial Ecology, QIIME 2 v2020.2, pipeline following the standard protocol for single-end libraries [44]. The operational taxonomic units (OTUs) table was generated using the DADA2 workflow [45]. Taxonomy classification of OTUs was performed with the SILVA 138 database [46] using the V-SEARCH algorithm. Output files generated by QIIME 2 pipeline were then visualized and analyzed within Phyloseq v1.30.0 [47], Ampvis2 v2.5.9 [48] and vegan v2.5-6 [49] R packages. Singletons, sequences unable to be classified at the kingdom level, and those classified as chloroplasts or mitochondria were removed from the data set before further analysis. All subsequent data analyses were performed at the genus level.
Alpha diversity measurements, including the number of observed OTUs, Simpson and Shannon indexes were estimated by rarefaction of 1324 sequences from each sample to adjust differences in library sizes across samples. The Kruskal-Wallis H test was performed to study significant differences in the alpha diversity indexes among groups. To measure the degree of uniqueness of each cover crop soil management concerning the variation of the composition of the community, the local contribution to beta diversity (LCBD) was calculated according to the procedure provided by [50]. Additionally, to explain variation in soil communities in terms of beta-diversity, a permutational multivariate analysis of variance (PERMANOVA) [51] with 999 permutations was calculated based on Bray-Curtis and Jaccard distance-based dissimilarity matrices. Differences among shared and unique OTUs at the different taxonomic levels among treatments were identified using proportional Venn diagrams generated with venn R package v1.9. To identify which taxa were differently abundant among the different treatments a differential abundance analysis as implemented in the DESeq2 R package [52] was performed.
Finally, a transformation-based redundancy discriminant analysis (tb-RDA) with the Monte Carlo permutation test based on 999 random permutations was used to determine which environmental variables (physicochemical, biological and functional properties of the soil) were the most significant in explaining the variation in bacterial community composition. The Hellinger data transformation was chosen to produce a more ecologically meaningful result, reducing the impact of low abundant species and double-zero problems [53]. To reduce the number of explanatory variables an automatic stepwise model building for constrained ordination methods was performed using the ordistep function with the forwarding selection implemented in the vegan R package. The tb-RDA model was tested for significant differences between the treatments using the anova.cca function and the significance of a relationship (p = 0.05) was assessed using the analysis of variance (ANOVA) based on 999 random permutations.

Erosion Measurements and Soil Physicochemical Properties
Runoff and soil loss were compared over an 8-year period in the plots managed with a gramineous cover crop (CC-GRA) and in the plot with conventional tillage (TILL). Data of the cover crop composed by a mix of different herbaceous species (CC-MIX) was only available for the last two years. For the whole period, in average terms, both cover crop managements showed effectiveness reducing cumulative soil losses (89%) and runoff (43%) ( Figure S2A,B) compared to TILL. For each year, runoff coefficients were moderately lower in the temporary cover crops (CC) treatments compared to TILL, ranging from 0.8-20.9% (6th and 8th years) and 1.4-24.6% (2nd and 1st years), respectively. Regarding soil loss, this trend was similar varying from 0.1 to 14.5 t·ha −1 (6th and 8th years) in CC and from 2.0 to 65.8 t·ha −1 (2nd and 5th years) in TILL. Only statistically significant differences among managements (p < 0.05) were found for annual runoff (6th year) and soil loss (4th and 5th years) ( Figure S1B).
Combined analysis of all soil physicochemical properties grouped the four treatments in three main clusters at a 69% cutoff value ( Figure 1A). The soil with a cover of natural vegetation (CC-NAT) for more than 20 years, which served as a benchmark, clearly clustered independently of the rest of the treatments. On the other hand, both treatments including a cover crop (CC-MIX and CC-GRA) were clearly different from the soil under tillage (TILL). In general terms, the analyzed physicochemical soil properties showed correlations between them (i.e., strong correlations were observed between pH and organic C, AL and CaCO 3 content or texture parameters and organic C) (Table S1). In all the cases, except for N content, soil properties differed significantly (p < 0.05) when soil managements were compared ( Figure 1B). Broadly speaking, whereas variables within the TILL treatment tended to show degraded soil characteristics, those of CC-NAT were the opposite, with intermediate values of those variables for CC-MIX and CC-GRA ( Figure 1B). Thus, the highest values for active limestone and CaCO 3 content and the lowest for organic C, exchangeable Ca, Mg and CEC were measured within the TILL management. This trend was the opposite at CC-NAT. An intermediate situation was found for CC-GRA and CC-MIX managements. Similar results were observed in most cases concerning these properties ( Figure 1B). Regarding organic N, no significant differences were found among managements, including CC-NAT where fertilization was null, although for this treatment the N values were the lowest ( Figure 1B).

Bacterial Diversity and Community Composition
After the quality control step of the NGS libraries, we obtained a total of 47,027 highquality sequences, with an average length of 300 bp. The DADA2 workflow produced 1245 unique OTUs, from which 1176 were identified as bacteria against the SILVA 138 reference database after removing unidentified (4.34%) and chloroplasts and mitochondria (1.20%) reads. A total of 1121 (95.3%) or 998 (84.9%) OTUs could be identified up to the genus or species level, respectively ( Figure S3). The rarefaction curves computed for each sample reached a clear asymptote for subsampling values above 1000, indicating that the samples had a good depth coverage, with samples under cover crop treatments saturating slightly slower than the TILL treatment ( Figure S4). Interestingly, all samples taken at the upper part of the slope showed much lower alpha diversity values for all indexes analyzed ( Figure S4).

Bacterial Diversity and Community Composition
After the quality control step of the NGS libraries, we obtained a total of 47,027 highquality sequences, with an average length of 300 bp. The DADA2 workflow produced 1245 unique OTUs, from which 1176 were identified as bacteria against the SILVA 138 reference database after removing unidentified (4.34%) and chloroplasts and mitochondria (1.20%) reads. A total of 1121 (95.3%) or 998 (84.9%) OTUs could be identified up to the genus or species level, respectively ( Figure S3). The rarefaction curves computed for each sample reached a clear asymptote for subsampling values above 1000, indicating that the samples had a good depth coverage, with samples under cover crop treatments saturating slightly slower than the TILL treatment ( Figure S4). Interestingly, all samples taken at the upper part of the slope showed much lower alpha diversity values for all indexes analyzed ( Figure S4).
All alpha-diversity indexes evaluated showed higher mean values in the soils with the presence of a cover crop (CC-NAT, CC-GRA and CC-MIX) compared to conventional tillage (TILL), with the soil with a natural cover crop showing the highest values, although those differences were not significant (p > 0.160) ( Figure S5). This effect was due to the All alpha-diversity indexes evaluated showed higher mean values in the soils with the presence of a cover crop (CC-NAT, CC-GRA and CC-MIX) compared to conventional tillage (TILL), with the soil with a natural cover crop showing the highest values, although those differences were not significant (p > 0.160) ( Figure S5). This effect was due to the differences obtained in the alpha diversity metrics across the terrain slope that were much lower in the upper part ( Figure S4).
In terms of beta diversity, bacterial community composition varied significantly among the different cover crop treatments ( Figure S6). Thus, PERMANOVA tests comparison, based on Bray-Curtis and Jaccard distances, showed that 39.48% (R 2 = 0.395; p = 0.016) and 35.32% (R 2 = 0.353; p = 0.006) of the total variance for both distance matrices, were significantly explained by the cover crop management in the orchard plots. The treatments with the presence of a cover crop of gramineous (CC-GRA) or spontaneous vegetation (CC-NAT) were differentiated apart from the CC-MIX and TILL treatments, which showed overlapping community composition ( Figure S6). In parallel, the LCBD values obtained for CC-NAT and CC-GRA were higher, in that order, independently of the taxonomic rank assessed (i.e., phylum, class, family or genus) (Figure 2).  Figures S7 and S8). A total of 590 bacterial species were identified; however, most of them belonged to uncultured or unidentified species. Consequently, data analyses were performed at the genus level. The core number of taxa detected simultaneously in all cover crops treatments were of 10 phyla, 17 classes, 35 orders, 45 families and 39 genera (Figures S7 and S8). As might be expected in soil ecosystems, the most abundant phyla, identified among all treatments, were Proteobacteria (36.15%), Actinobacteriota (27.89%), Firmicutes (16.52%) and Bacteroidota (14.99%) (Figure 2). From the 17 bacterial classes detected among all treatments, Alphaproteobacteria, Actinobacteria, Bacteroidia and Gammaproteobacteria were the most abundant, representing respectively 26.12%, 23.97%, 14.99% and 10.02% of the total (Figure 2, Figure S7). Five bacterial orders represented more than 50% of the total bacteria (Rhizobiales, Bacteroidales, Burkholderiales, Frankiales and Azospirillales, with seven bacterial orders (Dongiales, Elsterales, Lactobacillales, Christensenellales, Solibacterales, Gaiellales and Actinomycetales) being downregulated, showing significant differences among the different cover crops treatments due to their higher frequencies in the CC-GRA treatment ( Figure S9). At the family level, 45 families were detected across all soil samples, with Bacteroidaceae, Beijerinckiaceae, Geodermatophilaceae, Azospirillaceae, Xanthobacteraceae, Pseudonocardiaceae, Micromonosporaceae, Oxalobacteraceae, Lachnospiraceae and Lactobacillaceae being the most abundant (representing >55% of all bacteria; with mean frequencies ranging between 10.50% and 3.20%, in that order) (Figure 2, Figure S10). The relative abundance of nine families differed significantly among the different cover crop treatments, being the population for eight of those downregulated with higher population levels for the CC-GRA treatment, whereas for the Veillonellaceae family the trend was the opposite (Figure 2, Figure S10). At the genus level, a total of 39 genera comprised the core bacterial soil microbiome, being Bacteroides the most abundant (10.50%), followed by Microvirga, Blastococcus and Skermanella (>5.70%) (Figure 2). CC-GRA and CC-NAT were the cover crop treatments with the highest number of unique bacterial genera ( Figure S8). Twenty-eight bacterial genera differed significantly among the different cover crops treatments, with 15 genera being downregulated and 13 upregulated (Figure 2, Figure S11). Random forest classifier identified 10 bacterial genera as the most important, as being differentially expressed among the different cover crop treatments ( Figure S11). For example, Bradyrhizobium, Dongia and an unidentified Elsterales showed a higher log-relative normalized abundance in the CC-GRA treatment, whereas an unidentified Comamonadaceae, Lachnoclostridium and Microlunatus showed higher log-relative normalized abundance in the TILL and CC-MIX treatments ( Figure S11).

Soil Microbial Functional Activities
Soil functional indicators were greatly affected both by the presence of a cover crop and the plant species used. Combined analysis of enzymatic activities (API ZYM assay) or of the carbon-source utilization profiles (CLPP; Biolog EcoPlate assay) grouped the four treatments in two main clusters both for API ZYM ( Figure 3A) and the CLPP ( Figure 4A) assays. There was a clear trend to group the TILL and CC-GRA treatments together, whereas the treatment with a cover crop composed of a mixture of different plant species, either spontaneous (CC-NAT) or artificially sowed (CC-MIX), clustered together and independently of the former treatments ( Figures 3A and 4A). The functional indexes estimated for the API ZYM (TEA) and CLPP (AWCD Total and curves of increase in AWCD with time) were significantly higher (p < 0.045) for CC-NAT and CC-MIX than those obtained for TILL, with CC-GRA showing intermediate values. However, the number of enzymatic activities or carbon sources assimilated (S API ZYM , S CLPP ) was similar (p > 0.152) for all treatments ( Figure 3B and Figure 4B). In parallel, this higher functional activity of CC-NAT and CC-MIX was correlated with significantly (p < 0.041) higher numbers of culturable bacterial populations in CC-NAT and CC-MIX as compared to those measured in TILL, whereas CC-GRA showed intermediate values ( Figure 5). Soil respiration showed a similar trend than bacterial populations, with significant differences (p < 0.001) among all treatments, being highest for CC-NAT and lowest for TILL ( Figure 5).

Canonical Redundancy Analysis
A canonical transformation-based redundancy analysis (tb-RDA) was used to determine the existence of potential relationships between the bacterial community composition at the genus level and all the exploratory variables evaluated in the study including cover crop management, soil physicochemical properties and biological and functional soil activities. Overall, the first two tb-RDA axes (RDA1 and RDA2) explained 75.4% of the total variance in the composition of the bacterial community; with RDA axis 1 and RDA axis 2 explaining 46.6% and 28.8% of the variation in the bacterial community composition, respectively ( Figure 6). After forward selection, three variables including exchangeable Na (p = 0.007), pH (p = 0.011) and exchangeable Ca (p = 0.041) were selected as the best explaining the bacterial community structure. The RDA triplot ( Figure 6) shows that exchangeable Na (Na_exchange) plays an important role in the dispersion of the sites along the first axis, whereas exchangeable Ca (Ca_exchange) and pH do it in the second axis.  The position of some specific bacterial genera in the ordination indicates that the maximum abundance of most of those genera coincides with CC-NAT (e.g., Bacteroides, Massilia and Lactobacillus), CC-GRA (e.g., Bradyrhizobium, uncultured Acetobacteraceae and Xanthobacteraceae and Lactobacillus) and TILL (Geodermatophilus) soil managements ( Figure 5). This agrees with the previous results, indicating that those bacterial genera are favored by those specific cover crops (

Discussion
In this study we determined that the presence of a cover crop and its composition (regarding plant species) significantly influenced soil erosion and physicochemical quality indicators in an olive orchard under semiarid conditions. In parallel, bacterial community diversity and composition and the microbial functional activities in soil were greatly modified due to the presence of a cover crop, with a general trend to increase soil biological activities related to soil quality when compared to a tilled soil. However, these effects and their intensity differed greatly according to the plant species used. Figure S2B shows the annual relation between soil erosion and rainfall observing high runoff during the 1st, 7th and 8th years ( Figure S2A,B). This can be attributable to high annual precipitations during these years (858.4, 939.8 and 710.0 mm, respectively), greater than the mean of the rest of the years (489.9 mm) or the average annual precipitation in the area (534 mm). Average soil losses at the TILL treatment were 24.7 t ha −1 year −1 , being well above the maximum threshold of tolerable annual soil losses at Mediterranean catchments defined by several authors, e.g., 10 t·ha −1 , 5 t·ha −1 year −1 or 2 t·ha −1 year −1 [54][55][56], respectively. On the other hand, mean soil losses in the two treatments with the presence of a cover crop were much lower (2.7 t·ha −1 year −1 and 10.7 t·ha −1 year −1 , for CC-GRA and CC-MIX, respectively). Average runoff coefficients with CC-GRA and CC-MIX (6.3% and 19.9%, respectively) tended to be in the higher range of those measured in other studies in olive orchards at Southern Spain evaluating either spontaneous or gramineous cover crops [57,58]. Those studies reported average annual runoff coefficients of 1.2-6.8% and sediment yields of 0.8-1.2 t·ha −1 year −1 , with mean annual precipitations of 577-735 mm.

Effects of Cover Crops on Soil Erosion and Soil Physicochemical Quality Indicators
Our higher erosion records in the CC treatments could be explained mainly due to a combination of rainy years and a concentration of high intensity events in the period after the sowing in fall, in which the soil presented a very low ground cover. However, these results are in line with those measured in other studies in woody crops such as vineyards or almonds, under similar environmental conditions, where gramineous cover crops were implemented to control soil erosion processes, e.g., [59][60][61].
In general terms, our study has shown how temporary cover crops modified soil quality indicators and its beneficial effect on soil physicochemical properties compared to bare soil management based on tillage. Among the physicochemical parameters frequently used to characterize soil quality, soil organic carbon has been identified as one of the most discriminating indicators [62] together with cation exchange capacity, pH, available P, water storage characteristics and total N [63,64]. In our study, several of those physicochemical parameters ( Figure 1) were measured under the different managements, showing that in most of the cases all CC treatments have kept soil quality in an acceptable range of non-degradation, although being close to the lower limits proposed in [65] for organic olive orchards. Cover crops can be utilized in agricultural soils to manage C and N by increasing their content and altering their cycling and availability [66]. For the case of organic C, the TILL value was 0.61% (<0.84%, [65]) due to the herbaceous vegetation control in which soil was kept bare during the whole year. In our study, organic N in CC-MIX was 0.02% (<0.07%, [65]) despite the applied fertilization, showing high consumption by the cover crop species composing the mix. Nevertheless, cover crops in general terms, have a positive effect on soil organic carbon and nitrogen at woody crops contributing to the maintenance of soil quality but being strongly conditioned by specific crop species, agricultural management and climate conditions [67]. In summary, our results show beneficial effects on different soil physicochemical properties by cover crops but also that these effects take some years to be observed due to a combination of previous existent soil degradation and climate limitations [12,15,34,68]. Additionally, its impact can vary significantly from one farm to another [58,65], probably as a result of the differences at each site in the biomass production and species composition of the cover crop [69], as it has been shown as well for other woody crops in Southern Spain such as vineyards [68] or almonds [70]. All these results indicate the unsustainability of TILL and suggests that the use of temporary cover crops is an alternative treatment not only in terms of reducing soil erosion but also keeping soil quality regarding physicochemical properties, and favoring non-source pollution and nutrients losses, especially those associated with sediments [59].

Effects of Cover Crops on Soil Microbial Diversity and Composition
We found that the most abundant bacteria phyla present in the four managements were Proteobacteria, Actinobacteriota, Firmicutes and Bacteroidota. Although the proportions of those dominant phyla seem to keep relatively constant in different soils, the actual species present within a phylum may vary considerably among different soils and in terms of managements [71]. Our results indicated that the composition of the CC-NAT was different from the rest of treatments evaluated. Given the differences observed in the rest of soil variables, it is not surprising that we observe that the microbial composition of CC-NAT is slightly different from that of the other soil managements. Interestingly, phylum Firmicutes, mostly represented by the Bacilli and Clostridia classes, which include well-known antagonists and biocontrol agents, seem to thrive best in particular under the environmental soil conditions led by the presence of natural cover crops [72]. On the other hand, Proteobacteria (mainly represented by Rizobiales, Azospirillales and Burkholderiales) and Actinobacteriota are the dominant phyla in the rest of the sampled soils. The high abundance of Proteobacteria in soil is a good indicator of a high content of nutrients since most of them inhabit nutrient-rich soils [20,73]. In the same way, Actinobacteria is also considered an indicator of soil health [74], owing to members of this phyla playing important roles in the cycling of carbon, nitrogen, phosphorus, potassium and other elements. As soil saprophytes, they help to decompose organic matter into simpler molecules so they can be available to plants [75].
Some specific bacterial genera appeared as indicators of soil managements. Thus, Lactobacillus and Massilia were the most predominant genera for the CC-NAT soil. Lactobacillus spp. have been described as biological control agents against pathogenic bacteria and fungi [76] and plant growth promoting bacteria [77]. On the other hand, Massilia has been described in a wide range of habitats as a relevant component of the soil rhizosphere [78] exerting a relevant role in plant growth, mainly owing to the synthesis of secondary metabolites (i.e., auxin and siderophores) and by mobilization of elements such as phosphorus and heavy metals in phytoremediation processes [79], and also acting as a biocontrol agent against some soil-borne fungal pathogens [78]. Interestingly, Bradyrhizobium was the most abundant genus in CC-GRA despite no leguminous species being shown in this treatment and not being evident, while Bacteroids, Microvirga, Blastococcus and Skermanella were found in both CC-MIX and TILL, at similar proportions. Free-living Bradyrhizobium isolates lacking both nodulation and nitrogen fixation genes have been described in European soils [80], suggesting that this genus is extremely heterogeneous and has functions within the soil community that does not include symbiotic nitrogen fixation.
When analyzing the trends in soil bacteria diversity indexes, both in terms of richness (alpha diversity indexes, Figures S4 and S5) and on the composition and structure (beta diversity indexes, Figure 6, Figure S6), we found that those results were in agreement with the results obtained for soil erosion and physicochemical indicators. It is known that some agricultural practices that rely heavily on chemical fertilization, low plant diversity and chemical control of non-desirable weeds and pests adversely affect bacterial diversity [81,82]. As might be assumed from our results, the lowest bacterial diversity was expected to be observed in high disturbed soils (TILL), medium diversity in treatments with cover corps (CC-MIX and CC-GRA) and highest in low disturbed soils (CC-NAT). Cover crops provide an input of easily decomposable organic substrates through the plant residues and rhizodeposition available for the soil microorganisms [83,84]. These may explain the trend of finding higher values of alpha diversity indexes in soils with the presence of a cover crop, being highest for the CC-NAT ( Figure 6). Additionally, there was a trend towards finding greater alpha diversity in the lower slope zones (middle and lower) ( Figure S4), which may be explained by the accumulation of soil organic matter that is carried away from the surfaces of the higher zones by runoff.
Although, the CC-MIX treatment had very similar richness values to the CC-NAT treatment ( Figures S4 and S5), when analyzing the heterogeneity (beta diversity) of these bacterial communities, we found that these samples had a tendency to cluster with TILL rather than with the CC-NAT treatment, indicating that the composition and structure of CC-MIX soil biota was more similar to that found in TILL (Figures 2 and 6, Figure S6). In line with our results, [85] also reported that the main changes occurring on the soil bacterial community in Mediterranean olive orchard soil under different managements are basically linked to the bacterial community composition rather than to the presence/absence or number of soil microorganisms. The similar composition of the bacterial communities in TILL and CC-MIX treatments could be explained by the recent history of the CC-MIX plot, which has been under different soil treatment regimens (TILL and CC-GRA for six years and CC-MIX only for 2 years). Interestingly, although it seems that the mixed seed treatment is inducing an increase in bacterial richness, not enough time has elapsed for the bacterial communities to achieve the increased heterogeneity shown by CC-NAT. Nevertheless, this period seems to have been enough to induce relevant changes in the metabolic activities occurring on those soils as explained below.

Effects of Cover Crops on Soil Microbial Functional Activities
Our study revealed that the presence of cover crops (CC-MIX and CC-GRA) improved not only the physicochemical relevant properties such as organic C and cation exchange capacity but also several biological parameters related to soil quality (e.g., soil respiration, enzymatic activities and soil metabolic activity) compared to tillage management (TILL), reaching values closer to that of a physically undisturbed soil (CC-NAT) with a cover crop of herbaceous vegetation spontaneously present in the orchard.
Agricultural practices have been proven to modify soil microbiota allowing one to discriminate different soil managements in olive orchards. (e.g., [39,86]). However, the effect of the above ground vegetation and soil management in olive orchards on the bacterial population and soil microbial respiration is highly variable (e.g., [12,39]). Nevertheless, several authors (e.g., [26,39]) concluded that sustainable practices (e.g., use of cover crops, reduced and/or light tillage, pruning residues as a mulch, etc.) affect positively soil microbial diversity in Mediterranean olive orchards. In addition, there is a wide range of intermediate managements in which this answer is blurred. Several authors have asserted that bacterial community structure and diversity are strongly influenced by soil management, being evident that the positive effects of cover crops on bacterial biomass and enzymatic activity in Mediterranean areas [87,88]. In this regard, [24] compared several biological indicators under different soil managements in an olive orchard in South-Eastern Spain. Greater bacterial biomass and higher microbial functional diversity were found in ground covered soils compared to bare soils, highlighting the dominant effect of vegetation over other agricultural operations such as herbicide or mowing.
In Southern Spain, several studies [29,37,89] in organic olive orchards in two contrasting areas at the province of Córdoba have assessed the effect of different soil managements on bacterial population density, community composition and its activity. Those studies reported higher biological activity in grazing and mowing managements compared to tilled ones. In parallel, soil bacterial communities in all the managements showed a direct relation between aboveground (time of maintenance and vegetation diversity) and belowground bacterial diversity. Some authors [22,23] found higher enzyme activities in olive orchards under organic management with the presence of a cover crop as compared to paired conventionally managed orchards with chemical weed control. Nevertheless, site and sampling were determinants for detecting differences in these orchards located in South-Eastern Spain as well. In [90] the effects of cover crop management in an organic olive orchard on soil biological parameters at South-Eastern Spain were examined. The use of a cover crop (leguminous in this case) leaving residues on the soil surface drove to higher values of enzymatic activities, especially at the top 5 cm. It also seems clear in those studies that the cover crop residues used as a mulch and their subsequent decomposition result in higher values of chemical and biological parameters at these orchards under a Mediterranean climate.
In line with all the above-mentioned studies, we found a clear effect of the cover crop treatments on all the functional indicators measured including enzymatic and metabolic activities, soil respiration, bacterial population densities, with a trend to increase those variables as compared to the soil under tillage (TILL). Interestingly, the soil with the cover mixture (CC-MIX) reached values of those variables similar to those obtained in the soil that have maintained a cover of natural vegetation (CC-NAT) for more than 20 years, and which served as a benchmark.

Conclusions
This study describes the effects of four treatments based on the implementation of different ground covers (CC-NAT, CC-GRA and CC-MIX) and conventional tillage (TILL) on soil erosion, soil physicochemical and biological properties and soil microbial communities after 8 years of cover crop establishment. Our results showed that the use of a mixture of plant species (CC-MIX) or gramineous (CC-GRA) as a temporary cover in a non-organic olive orchard soil induced significant changes on soil indicators associated to soil quality when compared to a soil under tillage (TILL): i.e., reduced erosion processes, improved some physicochemical properties, modified the structure and diversity of soil bacterial communities and increased microbial functional activities. More importantly, some of the properties associated to poor soil quality can be reverted to a level similar to that of an undisturbed soil that had maintained a natural cover of spontaneous vegetation for decades (CC-NAT). All these results emphasize the need to promote the use of cover crops, in this case temporary, to enhance long-term soil conservation and sustainability of olive crop in Mediterranean sloping olive orchards.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/agronomy11071387/s1. Figure S1: Location of the trial and aerial view of the sampled managements (from Google Earth): TILL: tillage, CC-GRA: sown cover crop with gramineous, CC-MIX: sown cover crop with a mixture of species, and CC-NAT: cover crop with spontaneous vegetation. Figure S2: Erosion measurements during the 8-year period: precipitation and runoff (both in mm) and soil losses (t·ha −1 ) measured in the erosion plots. (A) Cumulative values and (B) Mean annual measurements with standard deviation. Different letters on the bars indicate statistically significant differences based in the ANOVA test (p < 0.05). Soil managements are: TILL: tillage, CC-GRA: sown cover crop with gramineous, CC-MIX: sown cover crop with a mixture of species. Figure S3: Taxonomic resolution of OTUs classified with the SILVA 13. 8 databases for bacteria. For each taxonomic classification, the lowest level was determined by plotting the total number of Reads and OTUs associated with the lowest level of classification. Figure S4: Comparison of rarefaction curves of observed OTUs obtained from an olive orchard soil under different cover crop managements (TILL: tillage, CC-GRA: sown cover crop with gramineous, CC-MIX: sown cover crop with a mixture of species) and CC-NAT (cover crop with spontaneous vegetation) and for each sampling point on the slope. The sampling points on the slope of the field for each treatment. Data shown are before rarefication. Figure S5: Alpha diversity comparison of an olive orchard soil under different cover crop managements (TILL: tillage, CC-GRA: sown cover crop with gramineous, CC-MIX: sown cover crop with a mixture of species) and CC-NAT (cover crop with spontaneous vegetation). The Kruskal-Wallis H test was performed to determine the existence of significant differences among the treatments. Data were rarified to 1324 sequences before analysis. Figure S6: Non-metric multidimensional scaling (NMDS) plot based on Bray-Curtis (left) and Jaccard (right) dissimilarity distances from an olive orchard soil under different cover crop managements (TILL: tillage, CC-GRA: sown cover crop with gramineous, CC-MIX: sown cover crop with a mixture of species) and CC-NAT (cover crop with spontaneous vegetation) for each sampling point on the slope. Stress values of the ordination analysis are shown at the top right of each plot. Figure S7: Prevalence Venn diagram showing the unique and shared bacterial taxa at different taxonomic levels obtained from an olive orchard soil under different cover crop managements (TILL: tillage, CC-GRA: sown cover crop with gramineous, CC-MIX: sown cover crop with a mixture of species) and CC-NAT (cover crop with spontaneous vegetation). The total number of taxa is indicated between brackets. * Number of taxa that stopped at a higher taxonomic level. Figure S8: Prevalence Venn diagram showing the unique and shared bacterial genera obtained from an olive orchard soil under different cover crop managements (TILL: tillage, CC-GRA: sown cover crop with gramineous, CC-MIX: sown cover crop with a mixture of species) and CC-NAT (cover crop with spontaneous vegetation). The total number of genera is indicated between brackets and the core and unique genera names for each treatment. * Number of OTUs that stopped at a higher taxonomic level. Figure S9: Differential OTU abundance analysis with DESeq2 R package at the order level. (A) Volcano plot representation of bacterial orders that are either up or down regulated; (B) Mean Decrease Accuracy plot of most important differentially expressed bacterial orders for microbial community obtained by random forest classifier; (C) Box-plot representation of the seven most abundant orders with a different significant abundance, for each taxon is represented the Log-relative normalised abundance among soil treatments and p-values from Wald test. Figure S10: Differential OTU abundance analysis with DESeq2 R package at the family level. (A) Volcano plot representation of bacterial families that are either up or down regulated; (B) Mean Decrease Accuracy plot of most important differentially expressed bacterial families for microbial community obtained by random forest classifier; (C) Boxplot representation of the seven most abundant families with a different significant abundance, for each taxon is represented the Log-relative normalised abundance among soil treatments and p-values from Wald test. Figure S11: Differential OTU abundance analysis with DESeq2 R package at the genus level. (A) Volcano plot representation of bacterial genera that are either up or down regulated; (B) Mean Decrease Accuracy plot of most important differentially expressed bacterial genera for microbial community obtained by random forest classifier; (C) Box-plot representation of the seven most abundant genera with a different significant abundance, for each taxon is represented the Log-relative normalised abundance among soil treatments and p-values from Wald test. Table S1: