Three-Year Investigation of Tillage Management on the Soil Physical Environment, Earthworm Populations and Crop Yields in Croatia

The aim of this study was to determine the environmental suitability of conservation tillage systems. A 3-year experiment was conducted in Croatia, to study the effects of different tillage treatments on soil properties, with the following: deep (DC), shallow tine cultivation (SC) and ploughing (P). Soil penetration resistance (SPR) was significantly greater in P compared to DC in all three years. In 2016, it was found at 30–40 and 40–50 cm; in 2017 at 10–20 cm; in 2018 at 0–10 and 10–20 cm. However, SC was significantly greater at 20–30, 30–40 and 40–50 cm compared to P and DC in 2017. The greater surface coverage in DC and SC (>30%) as compared to P (<1%) provided significantly higher soil moisture content (SMC) in maize (2016) and soybean (2018). In 2017, SMC in SC was significantly lower than in P and DC. Regarding all the 3 years, the agronomic structure in DC and SC had significantly greater crumb ratio compared to P, whereas P had significantly higher dust ratio than DC and SC. Throughout the 15 measurements, DC provided the most favorable soil habitat (11 occasions out of 15). In 2017, the earthworm abundance was significantly higher in DC compared to SC. In all the three years, DC resulted the highest yield, however the difference was not significant. Higher surface coverage and SMC positively impacted the ratio of agronomic structure (decreased dust and increased crumb ratio) and earthworm abundance. It can be concluded that DC and SC provided greater soil coverage which positively affected SPR, SMC, agronomic structure and earthworm abundance as compared to P.


Introduction
World population is gradually growing and reached 7.674 billion people in 2019 [1]. Such growth increases the demand for food, fiber and fuel. In this context, soil has been and will be one of the most valuable resources of mankind [2] since it is the support for 95% of food production [3]. To satisfy future population growth, it is necessary to preserve all productive soil, however, 33 million hectares of arable land is already severely degraded only in Europe [4]. Despite soil degradation, the negative effects of climate change on plant-soil relationship reduce food production [5,6].
Precipitation and air temperature, as the most important climatic factors, are under significant changes. According to climate change projections, precipitation decreases, and

Study Site and Description of the Experiments
The tillage experiment was set up on the border of Lukač (45°55′ N; 17°29′ E, 103 m a.s.l.) in the Pannonian region (Virovitica-Drava county) between 2015 and 2018 ( Figure  1). The area is mostly flat, with slope between 0° and 2° and land use cover was dominantly under croplands. The soil on experimental field belongs to Luvic Stagnosol (Siltic) with a silt loam texture [43] (IUSS Working Group WRB, 2015). According to depths, at depth 0-15 cm, soil is silty loam (28.8% sand, 49.5% silt and 21.7% clay), pH (H2O) 7.8 and organic carbon content is 2.23%, while at 15-30 cm depth, soil is silty loam (27.7% sand, 50.8% silt and 21.5% clay), pH (H2O) 8.0 and organic carbon 1.48%. During the 2015/2016 season, ploughing (P) and deep tine cultivation (DC) was studied, while in the next two seasons shallow tine cultivation (SC) was added as a treatment for evaluation. Each tillage treatment had three replicates. The individual plot area was 7500 m 2 (375 m long and 20 m wide). The tillage management data are shown in Table 1, while the data of crop production technology are shown in Table 2.  During the 2015/2016 season, ploughing (P) and deep tine cultivation (DC) was studied, while in the next two seasons shallow tine cultivation (SC) was added as a treatment for evaluation. Each tillage treatment had three replicates. The individual plot area was 7500 m 2 (375 m long and 20 m wide). The tillage management data are shown in Table 1, while the data of crop production technology are shown in Table 2.

Meteorological Data
The climate is continental with usually warmer autumn than spring. The average annual air temperature was 10.7 • C, and the average annual rainfall for the period 1965-1995 was 815.5 mm. During this time, the lowest measured average precipitation was 552.6 mm (1971), while the highest was 1114.8 mm (1972). Precipitation and temperature data (Virovitica Station) were obtained from the National Meteorological Service of Croatia ( Figure 2).   The study of the soil structure was performed according to the classification of Stefanovits (1992) [44]. Sampling was performed at the 0-10 cm depth in three replicates per treatment in each sampling period. The soil samples were air-dried and then carefully

Soil Sampling, Field Measurements and Laboratory Analyses
The soil conditions were measured five times per year over three growing seasons, from early spring to late autumn (Table 3). Soil penetration resistance (SPR) was measured with an electronic hand-pushed cone penetrometer (Penetrologger, Eijkelkamp, Netherland) using a cone with 1 cm 2 base area, with a 60 • included angle and 80 cm driving shaft. In total, 9 sampling points per treatment were measured at each sampling date and the data were grouped at 10 cm interval to a depth of 50 cm. Simultaneously to SPR measurements, the soil moisture content (SMC) was measured with a PT-I type measuring instrument (Kapacitív Kft., Budapest) in vicinity of each SPR measurement from 0-50 cm depth in 10 cm increment. The SMC was expressed in % (g g −1 ). The study of the soil structure was performed according to the classification of Stefanovits (1992) [44]. Sampling was performed at the 0-10 cm depth in three replicates per treatment in each sampling period. The soil samples were air-dried and then carefully sieved manually (60 shakes/min) on an agronomic sieve. Investigations of the ratio of different soil fractions such as clods (>10 mm), crumb (2.5-10 mm), fine crumb (0.25-2.5 mm) and dust (<0.25 mm) give the real state about the effect of tillage treatments on soil structure and the consequences of weather effects. The soil structure ratio can be determined from the weight of the fractions [18].
The soil surface covering was determined digitally with Adobe Photoshop CC 2019. All treatments were tested in nine replicates. The recordings were made after tillage as well as after sowing. The photos have a resolution of 2272 × 4608 and cover a minimum area of 6 m 2 . The soil surface covering results are shown in percentages.
The earthworm population was examined by hand-sorting in situ from a volume of 18.75 dm 3 (25 × 25 cm, 30 cm deep) in all plots in nine replicates [45]. The duration of hand-sorting lasted about 30-40 min, depending on the physical status of the soil. The examined sites were chosen randomly and the distance between soil blocks was about 5-10 m. The earthworms were sampled also like previous properties 15 times during three years ( Table 3). The earthworm abundance was expressed in ind m −2 [16]. Crops were harvested each year from each plot with harvester and the grain yield was weighed, and grain moisture was determined. Grain yields were presented with 0% moisture.

Statistical Analyses
Before the statistical tests were carried out, the datasets were checked for normality using Kolmogorov-Smirnov test. Soil coverage, earthworm abundance and crop yields follow the Gaussian distribution, while the SMC and SPR follow the normal distribution only after logarithmic transformation. Soil aggregates did not follow the normal distribution even after logarithm and Box-Cox transformation. A factorial ANOVA design was carried out to identify differences in the SPR and SMC (factors: tillage, sampling time and soil depth). One-way analysis of variance was used to examine the significant differences between the tillage methods and the soil cover, categories in the agronomic structure, crop yields and the number of earthworms. In case of a significant result of the analysis of variance, the groups with a significant difference were determined by Tukey HSD post hoc test at a p < 0.05. For soil aggregates, statistical comparisons were carried out with the non-parametric Kruskal-Wallis test (K-W). If significant differences at a p < 0.05 were observed, a multiple comparison of mean ranks post-hoc test was applied. Statistical data processing was performed using the IBM SPSS Statistics 25 software package (IBM, Armonk, NY, USA).

Influence of Tillage on Surface Cover
The soil cover data of the three-year experiment are presented in Figure 3. Soil cover is significantly lower in ploughing (P) compared to shallow tine (SC) and deep tine cultivation (DC) in each studied year. Surface cover in P is less than 1% of the soil each year. In DC surface cover was between 28.6% and 42.0%, while in SC it was between 36.1 and 46.0%.
processing was performed using the IBM SPSS Statistics 25 software package (IBM, Armonk, NY, USA).

Influence of Tillage on Surface Cover
The soil cover data of the three-year experiment are presented in Figure 3. Soil cover is significantly lower in ploughing (P) compared to shallow tine (SC) and deep tine cultivation (DC) in each studied year. Surface cover in P is less than 1% of the soil each year. In DC surface cover was between 28.6% and 42.0%, while in SC it was between 36.1 and 46.0%.

Soil Penetration Resistance
The soil penetration resistance (SPR) was significantly affected by tillage, year and depth and their interactions through all the three years of research (Table 4). Regarding the depths, SPR values were significantly higher with depths in all studied years. Depth × tillage interaction significantly increased SPR with depths at all treatments.
According to tillage, significantly higher SPR at P than at DC were observed in 2016 at 30

Soil Penetration Resistance
The soil penetration resistance (SPR) was significantly affected by tillage, year and depth and their interactions through all the three years of research (Table 4). Regarding the depths, SPR values were significantly higher with depths in all studied years. Depth × tillage interaction significantly increased SPR with depths at all treatments.   According to tillage, significantly higher SPR at P than at DC were observed in 2016 at 30-40 and 40-50 cm depth. In 2017, at 20-30, 30-40 and 40-50 cm depths, SC provided significantly higher values compared to DC and P. In 2018, P had significantly higher SPR values at all depths. At 0-10 cm and 10-20 cm depths, P showed significantly greater SPR than at DC. At 20-30 cm, P had significantly higher SPR than DC and SC. However, there were significant differences detected between P and SC at 30-40 and 40-50 cm.
According to measuring time, in 2016, significant differences between P and DC occurred in 4th and 5th measurements. In 2017, significantly higher SPR were observed at SC as compared to P and DC during 1st, 2nd, 3rd and 5th measurements. During 2018, significantly higher SPR were observed at P as compared to DC and SC at 3rd and 4th measurements. During the 5th measurement, P resulted significantly lower SPR than the other two treatments. Single factor analyses revealed significantly higher SPR at P than at DC in 2016 and 2018. In 2017, SC showed significantly higher SPR than other treatments, while during 2018 treatments significantly differed in the following order: P > DC > SC. Table 5 presents that all single factors as well as their interactions showed significant differences in soil moisture content (SMC) in all studied seasons.   Single factor analyses showed significantly higher SMC at DC as compared to P in 2016, while in 2017 and 2018 significantly higher SMC was recorded at DC compared to SC. In 2016, SMC was significantly greater (p < 0.001) in DC than P, at all depths and time of measuring (Table 5). In 2017, at 30-40 cm and 40-50 cm depth, SC recorded significantly lower SMC than the other treatments. During individual measurements, significant differences between treatments occurred in the 1st, 3rd, 4th and 5th measurements. In 2018, significantly higher SMC was detected at DC compared to P in all examined depths. Furthermore, DC showed significantly higher SMC compared to P (2nd, 4th, 5th) or SC (2nd, 3rd, 4th) in 2018.

Soil Structure
The average values of agronomic structure in tillage treatments for all years are shown in Figure 4.

Earthworm Abundance
Through the three years, the effect of tillage on earthworm abundance and its temporal dynamics in 15 in situ measurements were examined ( Figure 5). In 2016, higher number of earthworms were found in DC compared to P, but the difference was not significant (p < 0.05). According to temporal distribution, the highest earthworm abundance was counted in DC (72 ± 3 ind m −2 ) on 22nd July, while in P on The dust and clod ratios were the greatest in P in all the three experimental years, while the crumb ratio in P was significantly lower than in DC and SC. The fine crumb ratio was significantly lower in P in 2017 and 2018 as compared to DC and SC.
In 2017, the lowest ratio of clod fraction was measured in DC (13.12%), followed by SC (16.50%) and the highest (23.74%) in P, which significantly differed (p < 0.01) from the first two treatments. The significantly higher crumb ratio was determined in SC and DC as compared to P, while DC showed significantly higher fine crumbs as compared to DC.
In 2018, the highest clod (16.20%), dust (7.78%) and lowest crumb (37.38%) and fine crumb (38.64%) ratios were measured in P. Significant difference (p < 0.05) in clod, crumb and fine crumb fractions was found between SC and P, while in the dust fraction ratios showed significantly higher dust fraction at P as compared to DC (p < 0.05). Thus, in SC and DC significantly higher crumb (p < 0.01) and fine crumb ratio (p < 0.05) were recorded compared to P.

Earthworm Abundance
Through the three years, the effect of tillage on earthworm abundance and its temporal dynamics in 15 in situ measurements were examined ( Figure 5).
In 2016, higher number of earthworms were found in DC compared to P, but the difference was not significant (p < 0.05). According to temporal distribution, the highest earthworm abundance was counted in DC (72 ± 3 ind m −2 ) on 22nd July, while in P on 22nd June (48 ± 13 ind m −2 ). During the summer months (June, July and August) earthworm abundance decreased. Compared to August, in September more earthworms were counted in P (44 ± 8 ind m −2 ), and in DC (48 ± 23 ind m −2 ). This increase is due to an improvement in the quality of the habitat, i.e., temperature decreased, and there was 90.3 mm of rainfall in the meantime (Figure 2).

Earthworm Abundance
Through the three years, the effect of tillage on earthworm abundance and its temporal dynamics in 15 in situ measurements were examined ( Figure 5). In 2016, higher number of earthworms were found in DC compared to P, but the difference was not significant (p < 0.05). According to temporal distribution, the highest earthworm abundance was counted in DC (72 ± 3 ind m −2 ) on 22nd July, while in P on In 2017, in early spring (20th March) the highest earthworm abundance (28 ± 24 ind m −2 ) was found in P, followed by SC (16 ± 13 ind m −2 ), then in DC (12 ± 8 ind m −2 ). Regarding the temporal distribution, the highest earthworm abundance was recorded in May, more precisely in DC (152 ± 53 ind m −2 ), P (92 ± 20 ind m −2 ), then in SC (92 ± 33 ind m −2 ). In June, with the arrival of the warmer period, the earthworm abundance decreased in all treatments, i.e., −47.4% in DC, −34.8% in P and −47.8% in SC. This decrease resulted in a larger difference among treatments, thereby significant difference was found between DC and SC (p < 0.05).
In April 2018, more earthworms in P (24 ± 16 ind m −2 ) can be attributed to the faster warming of the topsoil. Through the three years, the highest number of earthworms was found in June, in DC (212 ± 59 ind m −2 ), P (156 ± 35 ind m −2 ) and SC (148 ± 24 ind m −2 ). The favorable change is due to the 148.8 mm precipitation between the two measurements and the good condition of the soybeans.
Regarding the temporal distribution, the largest earthworm populations were observed in June 2016, May 2017 and June 2018.

Crop Yield
The yields are presented in Mg ha −1 at dry moisture content. In 2016, the greatest maize yield was detected in DC (12.13 Mg ha −1 ), followed by P (11.91 Mg ha −1 ) ( Figure 6). In 2017, DC reached the highest yield (7.38 Mg ha −1 ), followed by SC (7.14 Mg ha −1 ) and P (6.90 Mg ha −1 ). After soybean harvest in 2018, there was no significant difference in yields between SC and P (4.20 Mg ha −1 ), while imperceptibly higher yield was obtained in DC (4.22 Mg ha −1 ). DC achieved the highest yields each year, while P reached the lowest or it was equal to SC.
The yields are presented in Mg ha −1 at dry moisture content. In 2016, the greatest maize yield was detected in DC (12.13 Mg ha −1 ), followed by P (11.91 Mg ha −1 ) ( Figure 6). In 2017, DC reached the highest yield (7.38 Mg ha −1 ), followed by SC (7.14 Mg ha −1 ) and P (6.90 Mg ha −1 ). After soybean harvest in 2018, there was no significant difference in yields between SC and P (4.20 Mg ha −1 ), while imperceptibly higher yield was obtained in DC (4.22 Mg ha −1 ). DC achieved the highest yields each year, while P reached the lowest or it was equal to SC.

Influence of Tillage on Surface Cover
The results of this study demonstrate how conservation tillage system (deep, DC and shallow tine cultivation, SC) can ensure proper soil cover. DC and SC showed significantly higher soil cover as compared to ploughing (P). Such findings are in agreement with other studies [13,21,46]. Surface cover is extremely important for soil physical quality, especially during the summer, as it plays a major role in soil and moisture protection [13].
Differences between surface cover between DC and SC can be largely attributed to tillage depth and soil mixing. This finding also corresponds to assessments by Raper's study (2002) [47], which reported that the depth of tillage also has a significant effect on stubble residue incorporation. This was also observed by Walther (2009) [48], with a higher ratio of soil surface remaining covered by shallower tine cultivation.

Soil Penetration Resistance
In 2016, soil penetration resistance (SPR) was statistically higher in P than DC at 30-40 and 40-50 cm. In 2017, significantly higher SPR was obtained in SC at 20-30, 30-40 and 40-50 cm compared to P and DC. Tillage and their working depth had a variable effect on SPR. Bogunovic et al. (2018) [49] and Ren et al. (2019) [50] observed similar research findings. In addition, SPR values increased with depth to more than 4 MPa, especially from about 30 cm and below. Birkás et al. (2018c) [51] found the occurrence of the compact layer (plough pan) at~28 cm depth. The SPR in P was significantly higher than DC (0-10 and 10-20 cm) in 2018 in this study.
Residue retention under DC contributed to water conservation through multiple effects on water regime which also has effects on SPR. Very low percentage of stubble residues in P [16,18] may negatively affect soil physical properties in several ways, such as higher soil temperature [52], higher soil evaporation [53], lower soil moisture content (SMC) [54] and greater deterioration of soil structure and soil settling [55] through the kinetic energy of raindrops. However, Afzal et al. (2008) [56] found that rapid increase of soil temperature in bare soil, especially in early spring, might be recognized as the crucial factor for faster seed germination.
According to the measurement time, in 2016 and 2017, the SPR values in the 1st measuring time were significantly lower than all the others. Despite the 170.5 (2016) and 148.9 mm (2017) precipitation between the 1st and 2nd measurements, the 2nd measurement time resulted significantly greater SPR. An additional 47.7 (2016) and 48.0 mm (2017) precipitation from the 3rd to the 4th measurements contributed to a significant difference in SPR. There was significant difference between 4th and 5th measurements in 2017, due to the increase of temperature and low precipitation (29.8 mm). In 2018, the significantly lowest SPR value was in the 2nd measurement which was due to tillage operation, and it significantly differed from all the other measurements.
Considering the measurement time, the impacts of climatic conditions (heavy rainfall and its occurrence), but also water consumption are key factors. These changes in the soil are best demonstrated by experiments with artificial irrigation [57]. They highlighted that on the top of positively affecting productivity, it encourages a progressive modifying physical parameter, emphasizing greater compaction and humidity, and with lower porosity and soil weight.

Soil Moisture Content
In the present work conservation tillage treatments (DC and SC) resulted in relatively higher SMC in all years and measurement time as compared to P, which supported the hypothesis in this study. The SMC values of P treatments in 2016 at 40-50 cm depth were the lowest (12.4%). The SMC values drastically decreased in the two lower layers, which is closely associated with plough pan. In 2017, the SC as the shallowest tillage treatment gave significantly lower SMC compared to DC and P at 30-40 and 40-50 cm depth. In 2018, the SMC values in all depths were significantly greater in SC than in P. Furthermore, SMC values did not fluctuate so strongly in DC and SC as compared to P. Our study indicates that the ratio of surface cover has a positive relationship with SMC in DC and SC treatments, especially at the surface layers. Similar findings were also reported by Kalmár et al. (2013) and Dekemati et al. (2019) [13,16]. Dekemati et al. (2019) [16] found that at 10 cm depth, a significantly higher and the greatest SMC were measured in DC (26.06%) as compared to P (21.24%) which was the driest.
Current agricultural production in the Central European region is faced with unbalanced precipitation, therefore the water availability in soil is a key factor in dry farming which influences the growth and development of crops. In areas where irrigation is not possible, soil conservation practices, such as minimum disturbance and higher surface residue are highly valued as soil moisture management tools. Many papers had been published about the favorable effects of reduced disturbance and higher surface residue with the aim of better control and moisture conservation [14][15][16]58].

Soil Structure
Nowadays, due to frequent heavy rains and prolonged droughts, the protection of the soil surface with plant cover has become important. In our experiment, the highest ratio of clod was reached by P which is in agreement with a similar study in Hungary [59]. Weather, i.e., too high or too low temperature, strong rainfall accompanied by a strong wind and ploughing as a tillage method provide unfavorable circumstances for soil health, respectively. In addition to this, many papers had been published about positive effects of mulch residues on the soil [11,12,18,20]. The present results confirm previous works as DC provides the lowest ratio of dust. With similar results, Baker et al. (2005) [35] drew attention to the fact that conservation tillage achieves a 50-90% reduction in dust fraction, especially when the number of tillage interventions is reduced. However, crumb fraction ratios between DC and SC treatments showed no difference, while both treatments showed significantly higher crumb ratio than P. This result is not only due to the surface cover which is provided by DC and SC, but also because the other favorable soil characteristics of these tillage interventions, e.g., increased SMC, greater earthworm abundance, more intensive plant root growth. In contrast to P, moderate disturbance at SC and DC, and the incorporation of plant residues into the topsoil result in favorable dynamic biological activity. A similar phenomenon was observed by Kalmár et al. (2013) and Bogunovic et al. (2019) [13,19] in their relevant experiments.

Earthworm Abundance
This study showed that tillage did not have significant impact on earthworm abundance, except for one case (out of 15) on 21 June 2017 when DC had significantly higher abundance compared to SC. Considering the hypothesis of this work, conservation tillage significantly changed the soil conditions (greater SMC, higher ratio of crumb and lower ratio of clod and dust) which had some impacts on abundance. Which means that the hypothesis was partially justified. The earthworm abundance was greater in 11 occasions out of 15 in DC. Similar results were found Bogužas et al. (2010) [60], who compared shallow cultivated soil to P, and found 51% more earthworms in the shallow cultivated soil. In this study, at the two early spring dates (20 March 2017 and 24 April 2018), the earthworm abundance was greater in P, very likely due to higher soil temperatures in this treatment. In P, the higher earthworm population is presumably attributable to the uncovered surface, which warms up faster in the spring. In line with this, Eriksen-Hamel and Whalen (2006) [61] found that abiotic factors, such as temperature and SMC, largely regulate the dynamics of earthworm populations.

Crop Yield
This study showed that conservation agriculture (CA) treatments (especially DC) did not result in significantly higher yields. Among the three treatments, in all the three crop years, P reached the lowest yields. Similar results were obtained by Bescansa et al. (2006) [62], who pointed out that tillage or crop residue did not have significant impact on yields. However, in dry years, ploughed treatment resulted in lower yield. In line with this, Bogunovic et al. (2020) [24] reported that crop yields increased under loosened soil compared to ploughed. Slightly higher yields in DC could be related to better retention of water through the observed soil physical properties.

Conclusions
Mulch covered treatments (deep, DC and shallow tine cultivation, SC) provide higher topsoil protection and favorable soil physical circumstances indicated by generally lower compaction, and higher water conservation as compared to uncovered ploughed soil according to this study. The soils treated with deep tine cultivator had the lowest values of soil penetration resistance (SPR), while ploughing (P) created a compacted soil environment for plant growth. Ploughing retained less water than DC. The outcomes of this study highlighted that the use conservation tillage method (DC and SC) is more beneficial than conventional ploughing, since it increases residue cover and proportion of macro-and mezo-aggregates. The largest earthworm populations were observed in DC (11 occasions out of 15). The most favorable edaphic conditions for earthworms were provided in DC in June 2016, May 2017 and June 2018. Based on these results, it can be concluded that ploughing under the studied agro-ecological conditions should be avoided since it increases the proportion of micro-aggregates and dust, and decreases the crumb ratio. It can be assumed that changes in soil physical and biological conditions are much faster, however, the effect of these changes on the yield improvement manifests slower. Considering the erratic weather conditions, the storing of water in the soil is extremely important which can be achieved by appropriately chosen soil tillage system (in this case it was represented by DC and SC).