Is Mechanized Harvesting of Shrubs Grown for Energy Purposes Environmentally Sustainable?

: Mechanized harvesting of shrub formations as part of sustainable forest management not only helps reduce the risk of forest ﬁres in Mediterranean environments but also provides economic beneﬁt from the extracted biomass, contributing to the development of the bioeconomy. However, these mechanized operations require an environmental impact assessment to identify the short-term impacts, both positive and negative. This is especially important in the Mediterranean basin given the speciﬁc climatic conditions which exist. In this study, the environmental impact of mechanized shrub harvesting is analyzed in relation to (i) changes in biodiversity and in the presence and growth of species; (ii) physical and chemical properties of the soil; and (iii) changes in forest ﬁre risk. For this purpose, a pre-harvest inventory was conducted and post-harvest monitoring schedules of 1- and 2-year durations were established in three characteristic Mediterranean shrubland formations located in the northern–central area of the Iberian Peninsula. Our results reveal that the recovery rates in biodiversity indices after harvesting were very high, with values ranging from 30 to 70% depending on the site. Two years after harvesting, the species coverage was similar to the pre-harvest scenario in some locations, although not with regards to height, the ericaceous species being those with the greatest sprouting capacity. Signiﬁcant changes in the physical and chemical properties of soils were also observed. In this regard, negative impacts such as soil compaction or slight acidiﬁcation were identiﬁed at some sites. However, positive effects were also found such as an increment in carbon and nitrogen content after harvesting, along with increased litter quantity a year from the clearing operation due to biomass residue left on the ground after harvesting. Furthermore, mechanical harvesting effectively modiﬁed ﬁre behavior in all the shrub formations 2 years after clearing, with a notable reduction in ﬁre risk at all the studied sites.


Introduction
Shrubland comprises dense thickets of evergreen sclerophyll shrubs and small trees. It is a biome that extends across the temperate zones of the Earth with a Mediterranean climate, characterized by low precipitation and a marked dry season [1]. In particular, these shrublands are found in the mid-latitudes of California, southwest and southeast Australia, Chile, the Western Cape of South Africa, and the Mediterranean basin. The distribution and composition of shrubland in these regions is determined by the combined influence of climate, soil substrate, water availability, and disturbance regimes [2], with both successional and climatic shrub formations being differentiated according to the greatest determining features in their composition. Shrubs and bushes are highly abundant in the Mediterranean basin, accounting for more than 29% of the main forest habitat types [3], mainly concentrated in Portugal, Spain, France, Italy, and Greece. They take different forms depending on geographical location, altitude, exposure, and soil composition, and are mainly associated with with a temperate climate and a prolonged summer drought period. In this regard, it is important to determine the temporal changes in biodiversity and species dominance, as well as the effect on soil properties, erosion, and forest fire risk, since these aspects will differ greatly from other regions due to the specific climatic conditions of the Mediterranean basin. The majority of studies have only focused on tree vegetation [17] and the approaches developed to assess short-term impacts are unsuitable for addressing the complex, rapid changes which occur in a Mediterranean shrub ecosystem. A clear understanding of all aspects involved will also allow an optimal rotation period to be established for each shrub formation. Biomass harvesting can then be scheduled, and the economic profitability analyzed. Thus, the productivity and sustainability of the shrubland will not be compromised.
The objective of this study was to evaluate the main environmental effects of mechanical harvesting in shrubland communities under Mediterranean conditions. The specific objectives were to test the impacts of harvesting on (i) biodiversity and regeneration capacity of the species; (ii) certain physical and chemical properties of soils; and (iii) certain factors linked to forest fire risk. The working hypothesis was based on considering all of these impacts together over a short (one-year and two-year) time period after the mechanical harvesting.

Study Area and Selection of Sampling Points Prior to Shrub Harvesting
The evaluation of harvesting impacts was carried out at three locations situated in the northwest of the Iberian Peninsula (Table 1). The selected shrub formations were characteristic of Mediterranean shrubland. The machines used for harvesting operations were two different commercial harvesting systems: (i) one baling harvester (BIOBALER WB55; prototype) which produced bales and (ii) one mulching harvester (RETRABIO, own design), that produced crushed biomass. In each of the study areas where harvesting was planned, we identified locations with a medium shrub coverage, using aerial photographs. Edge effects were avoided. A 50 × 50 m grid was then overlaid on each of the selected locations using a GIS tool. Suitable sampling points for performing the environmental impact assessment were determined on the grid (six per location). The climatic, edaphic, and physiographic characteristics of the study sites are shown in Table 1. The main shrub species in terms of abundance were Genista cinerascens Lange, Genista florida L., Halimium lasianthum (Lam.), and Erica australis L. Other shrub species, although less abundant, were also represented ( Table 2).  In the field, with the aid of the Global Positioning System (GPS), a permanent plot was established at each of the six sample points identified per location (using real coordinates). Based on the reference point coincident with the coordinates, a base line of 18 m in length was established at each of the sample points, along which transects were established to inventory vegetation and perform soil and litter sampling.

Vegetation Transects and Soil and Litter Sampling
For the purpose of the vegetation inventory, three 25 m long transects perpendicular to the base line were established in each plot. The number of inventory points per location was 18 in the case of vegetation and 36 for soil and litter sampling. However due to the impossibility of clearing some of the plots as planned at the beginning of the activity (pre-harvest scenario), the number of inventories was reduced post-harvesting. Table 3 shows the final inventories at each stage and location. In each transect, all species were identified, the beginning and end of each individual plant crown were determined following the line intercept method [21], and their heights were recorded. For non-woody plants, we differentiated between annual and perennial grasses, as well as the presence of mosses.
Additionally, a point 6 m along the (18 m) baseline was used as a reference point for litter and soil sampling. For the purposes of litter sampling, square structures of 50 × 50 cm and 25 × 25 cm were placed 1 m either side of the reference point. We differentiated the litter layer and recognizable remains (L), which were collected within the larger square, and the fragmented and humified remains (FH), which were collected inside the smaller one. We also noted the depth of each stratum. All these samples were then placed in separate plastic bags for subsequent transfer to the laboratory. Within each of the square structures, a metal cylinder of 6.5 cm in diameter was then inserted into the superficial (0-5 cm), subsuperficial (5-10 cm), and deep (15-20 cm) horizons to perform the soil sampling. The number of soil samplings and litter sampling at each stage of the environmental assessment are shown in Table 3. These samples were then used for the analysis of the different variables studied.

Coverage, Height, and Biodiversity Analysis
In a first step, an analysis of the total coverage was performed using the pool of the main woody species at each location. Coverage was calculated in linear meters and expressed as a percentage. These percentages were determined for each transect of each site, and the average value obtained was used. We then compared the changes in the percentage cover of the pool of species over time and for each location using a two-way ANOVA with an unbalanced design due to the differences in sampling points before and after harvesting. To normalize the error, percentage cover data was arcsine-transformed. The significance of the means was tested using the Newman-Keuls test. A second analysis of the coverage was then performed separately at the woody species level for the main species at each location using the same statistical criteria as those for the total coverage. The heights at species level were also analyzed. The effect of time elapsed since harvesting on height was also analyzed using a two-way ANOVA with an unbalanced design.
Indices of richness (Margalef), diversity (Shannon), and dominance (Simpson) were also calculated to evaluate the effect of shrub harvesting on the biodiversity at each location. We considered all the woody species separately, whereas the herbaceous species were considered as a whole. The Margalef index [22] of ecosystem richness only considers the total number of species found at each sampling, including species belonging to the regeneration stratum. In addition to the total number of species, the Shannon index [23] estimates the diversity of the ecosystem, taking into account the proportion of individuals of a given species versus the total individuals present in the community analyzed. Finally, the Simpson index [24] indicates the degree of dominance present in the ecosystem, taking a certain number of species present in the habitat and their corresponding relative abundance.
These indices were calculated independently for each sampling point and transect, and their average value was used to calculate the index of each location. The following mathematical expressions were used: where S is the number of species present and N is the total number of individuals in the sample; where p i is the proportion between the individuals of one particular species (i) and the total number of individuals found (N); and where n i is the number of individuals found of one particular species, and N is the total number of individuals found of all the species. For both the Margalef and Shannon indices, values below 2 are considered areas of low diversity, values between 2 and 3 indicate average biodiversity, and values greater than 3 are indicative of high biodiversity [22,23]. Simpson's dominance index provides results between 0 and 1, where values close to 1 indicate the dominance of one species over the others, i.e., more homogeneous ecosystems [23].

Soil Variables
The following variables were calculated prior to harvesting and again one year and two years after clearing: (i) bulk density (D a ; g cm −3 ), (ii) total carbon content (C; %), (iii) total nitrogen content (N; g kg −1 ), and (iv) pH (KCl).
To determine soil bulk density, one undisturbed soil sample (from a steel cylinder of 5 cm height and 6.5 cm diameter) was taken for each of the three sampling depths at every sampling location. The remaining soil samples were then air-dried, weighed, and sieved (<2 mm) to obtain the fine soil proportion of each sample, on which the chemical analyses were performed, as well as the particle size distribution to obtain the soil texture, which was determined using the pipette method (UNEP-UN/EC Method 9100SA). The pH was potentiometrically measured (Metrohm mod 691) in the supernatant suspension of the sieved soil samples in 1N KCl (soil/solution ratio 1:2.5) after an equilibration time of 30 min. The total carbon (C) and nitrogen (N) contents of the grounded soil samples were analyzed using a dry combustion method with a LECO elemental analyzer (model CN-2000, LECO Corporation, Saint Joseph, MI, USA). C and N concentrations were calculated at 65 • C.
A textural analysis of each horizon was carried out at some locations as a consequence of the results obtained with regard to carbon and nitrogen contents. This analysis was performed using the pipette method. Additionally, and as a consequence of the results obtained for the pH, a cation exchange was also performed. The concentrations of the exchangeable basic cations sodium, potassium, calcium, and magnesium and the exchangeable acid cations iron, manganese, and aluminium were determined in a 0.1 M barium chloride extract of the soil. The soil extracts were obtained by shaking 2.5 g air-dried soil (<2 mm) with 30 mL barium chloride solution for 2 h, centrifuging at 3000 rpm for 10 min, and vacuum filtering the supernatant liquid (0.45 lm membrane filter) (reference methods ISO 11260 and ISO 14254). The exchangeable cations in the extract were determined by ICP-OES (Perkin-Elmer, Krakow, Poland, Optima 5300 DV).
The differences between years (before and after harvesting) among the soil variables were analyzed using the nonparametric Kruskall-Wallis test. Mann-Whitney U tests (non-normality variables) were carried out to determine whether significant differences existed between repetitions over time in the cation exchange analysis. Additionally, Pearson's correlation coefficient (r) and the significance level (p) were used to determine the level of similarity among certain variables. All the statistical analyses were performed using the statistical package IBM SPSS STATISTICS 20.0 (IBM Corp., SPSS Statistics for Windows, Armonk, NY, USA), and the level of significance was set to p = 0.05 for all the analyses.

Evaluation of Fire Risk through Simulation Modelling
To characterize and quantify changes in the combustibility of shrub formations subjected to mechanical clearing, custom fuel models were simulated for each of the locations at three separate moments in time: before harvesting, a year after harvesting, and two years after harvesting. For this purpose, three whole, representative plants of the main shrub species at each site and date of inventory were selected, as well as the biomass contained in 3 replicates of 1 m 2 each at each location. The maximum dominant height of the main species of the collected samples was also recorded.
Once in the laboratory, the samples were dried at 65 • C to constant weight and fractionated, differentiating between live material <0.6 cm (shoots and leaves) and dead material of the three defined size classes (<0.6, 0.6-2.5, and 2.5-7.5 cm).
Fire behavior simulations were performed using a software package specifically designed for this purpose (Behave Plus 5.0). We considered a different simulation scenario for each study site and data (before and after harvesting). Four variables of interest were calculated: (i) fire propagation speed (S r ; m min −1 ); (ii) heat per unit area (H a , kg m −2 ); (iii) fire line intensity (F i ; Kw m −1 ); and (iv) flame length (F l ; m). Simulations were carried out considering three possible wind scenarios (mild, moderate, and strong), as well as considering the average slope value and the physiographic and climatic characteristics of the studied areas.

Impacts of Harvesting on Coverage, Height, and Species Biodiversity
The analysis of total coverage ( Figure 1) revealed a significant reduction during the first year after harvesting at all study locations. However, two years after the clearing operations, locations FA and FB showed no significant differences compared to the pre-harvest scenario (Figure 1), while a significant reduction in total coverage was still evident at location NV (Figure 1).

Impacts of Harvesting on Coverage, Height, and Species Biodiversity
The analysis of total coverage ( Figure 1) revealed a significant reduction during the first year after harvesting at all study locations. However, two years after the clearing operations, locations FA and FB showed no significant differences compared to the pre-harvest scenario (Figure 1), while a significant reduction in total coverage was still evident at location NV (Figure 1). The analysis of coverage at the level of the most important woody species at each location ( Figure 2) revealed that most of the species showed significant reductions in their coverage both one year and two years after harvesting. However, this was not the case for ericaceous species, i.e., Erica australis and Calluna vulgaris. These showed significant reductions one year after harvesting but also significant increases in their coverage two years after harvesting compared with the pre-harvest scenario ( Figure 2). The analysis of heights of the shrub species ( Figure 2) before and after harvesting revealed that in all cases (including ericaceous species) there were significant reductions both one and two years after clearing compared with the pre-harvest scenario.  The analysis of coverage at the level of the most important woody species at each location ( Figure 2) revealed that most of the species showed significant reductions in their coverage both one year and two years after harvesting. However, this was not the case for ericaceous species, i.e., Erica australis and Calluna vulgaris. These showed significant reductions one year after harvesting but also significant increases in their coverage two years after harvesting compared with the pre-harvest scenario ( Figure 2). The analysis of heights of the shrub species ( Figure 2) before and after harvesting revealed that in all cases (including ericaceous species) there were significant reductions both one and two years after clearing compared with the pre-harvest scenario.

Impacts of Harvesting on Coverage, Height, and Species Biodiversity
The analysis of total coverage (Figure 1) revealed a significant reduction during the first year after harvesting at all study locations. However, two years after the clearing operations, locations FA and FB showed no significant differences compared to the pre-harvest scenario (Figure 1), while a significant reduction in total coverage was still evident at location NV (Figure 1). The analysis of coverage at the level of the most important woody species at each location ( Figure 2) revealed that most of the species showed significant reductions in their coverage both one year and two years after harvesting. However, this was not the case for ericaceous species, i.e., Erica australis and Calluna vulgaris. These showed significant reductions one year after harvesting but also significant increases in their coverage two years after harvesting compared with the pre-harvest scenario ( Figure 2). The analysis of heights of the shrub species (Figure 2) before and after harvesting revealed that in all cases (including ericaceous species) there were significant reductions both one and two years after clearing compared with the pre-harvest scenario.  Species diversity (Shannon index), species richness (Margalef index), and species dominance (Simpson index) for each of the study locations prior to harvesting, as well as one and two years after harvesting, are shown in Table 4. The general pattern observed one year after shrub harvesting at all locations was a decrease in the number of species, which resulted in reductions in the Margalef and Shannon indices. These reductions were particularly noteworthy at site FB (42% and 67% depletion respectively in these indices; Table 4). The reduction in species was least noticeable at NV (reductions of around 24% and 13% in relation to the Margalef and Shannon indices, respectively). FA showed an intermediate behavior compared to the other two sites, with a decrease in the Shannon index of 37%, but a more pronounced decrease in the Margalef index, i.e., 60% (Table 4). However, two years after shrub harvesting, the reductions in the Margalef and Shannon indices compared to the situation pre-clearing were far less pronounced at all the sites (8-11% for NV, 3-14% in FB, and 6-33% for FA, Table 4). Furthermore, in some cases the observed values were even closer to the pre-harvest value, or slightly higher in the case of FB (Table 4).
With respect to the Simpson Index (Table 4), there was a significant increase in values at all the studied sites during the first year after harvesting, indicating that the shrub communities were less diverse and that the degree of dominance had increased. This increase was particularly noteworthy in FA (increments of 79%), followed by FB (68% increment) and NV (58% increment). However, the values for this index were lower for all sites in the second year after harvesting, reaching the pre-harvest values in the case of FB, although such values had not been reached at NV or FA.

Impacts of Harvesting on Physical and Chemical Properties of the Soil
The average values for the main soil variables studied before and after shrub harvesting for each of the three soil horizons analyzed are shown in Table 5.
Shrub harvesting at site NV led to a significant reduction in soil pH one year after clearing, both in the superficial and subsuperficial horizons (p < 0.05), as well as in the deep horizon (p < 0.01). However, after two years, the soil pH reached very similar values to those registered before clearing, with no significant differences (Table 5). In addition, there was a significant increase in the bulk density of the soil both one year and two years after harvesting, although only in the superficial horizon (p < 0.05).
Mechanical harvesting at site FA led to a significant reduction in the nutrient content (N, C) in the superficial horizon (p < 0.05) one year after harvesting (p < 0.05), while in the other two horizons the opposite effect was observed, i.e., a significant increase in C and N, which was especially significant in the deep horizon (p < 0.01). However, two years after harvesting there has been a significant increment in the nutrient content (N, C) in all three horizons (p < 0.01). A significant increase in soil bulk density (p < 0.01) was also observed at this site one year after harvesting along with reduced pH (p < 0.01) in the deep horizon (Table 5).
Harvesting at site FB led to a significant increase in C and N (p < 0.05) for both years, as well as a decrease in pH (p < 0.05) in the deep horizon one year after harvesting, with no changes in the two upper horizons.  (2015), as well as one year (2016) and two years (2017) after clearing: bulk density (D a ; g cm −3 ), nitrogen (N; g kg −1 ) and carbon (C; %) content, pH (KCl).  With respect to the litter analysis (Figure 3), all the sites showed significantly higher values for the litter layer and recognizable remains (L) one year after the mechanical harvesting compared to the pre-harvest scenario. However, no significant changes were identified for any of the sites in terms of fragmented and humified remains (FH) between the pre-and post-harvesting scenarios (Figure 3).
With respect to the total litter thickness, the three studied sites showed a significant increase in this variable one year after harvesting. In contrast, while no significant differences were found between the total litter thickness two years after harvesting and the pre-harvest situation at NV and FB, at FA it still was significantly greater (Figure 3).
With respect to the litter analysis (Figure 3), all the sites showed significantly higher values for the litter layer and recognizable remains (L) one year after the mechanical harvesting compared to the pre-harvest scenario. However, no significant changes were identified for any of the sites in terms of fragmented and humified remains (FH) between the pre-and post-harvesting scenarios (Figure 3).
With respect to the total litter thickness, the three studied sites showed a significant increase in this variable one year after harvesting. In contrast, while no significant differences were found between the total litter thickness two years after harvesting and the pre-harvest situation at NV and FB, at FA it still was significantly greater (Figure 3). The textural analysis for each horizon revealed that location FA showed a significant increase in clay, with values ranging from 18% to 30% depending on the horizon, while these values were even higher at location FB (25-40%).
The cation exchange capacity analysis (Table 6) showed a significant increase in aluminium cations at all the sites and for all soil horizons, in which a significant acidification was observed during the first year after harvesting. However, no significant changes were found after two years. Moreover, a significant increase in calcium was observed in the deeper horizons (15-20 cm) at sites FB and FA, where a significant reduction in pH was also recorded during the first year after clearing. However, no significant changes were found for the rest of the basic cations (Na, K, Ca, Mg), nor the exchangeable acid cations Fe and Mn (data not shown), either for the first or the second year post-harvest.

Impacts of Harvesting on Fuel Load and Fire Risk Reduction
The pre-and post-harvest fire parameters for each of the sites are shown in Table 7. Before harvesting, site FA showed higher values for all the analyzed variables, indicating increased fire risk compared with NV or FB. This is due to the high average height of the broom plants in this area, as well as the presence of an abundant herbaceous load, resulting in greater fire propagation speed, heat per unit area, fire line intensity, and flame length than at the other sites ( Table 7). The effectiveness of shrub harvesting in reducing the risk of fire was greatest at site NV. Thus, considering an intermediate scenario with a wind speed of 15 km h −1 , the fire propagation speed was The textural analysis for each horizon revealed that location FA showed a significant increase in clay, with values ranging from 18% to 30% depending on the horizon, while these values were even higher at location FB (25-40%).
The cation exchange capacity analysis (Table 6) showed a significant increase in aluminium cations at all the sites and for all soil horizons, in which a significant acidification was observed during the first year after harvesting. However, no significant changes were found after two years. Moreover, a significant increase in calcium was observed in the deeper horizons (15-20 cm) at sites FB and FA, where a significant reduction in pH was also recorded during the first year after clearing. However, no significant changes were found for the rest of the basic cations (Na, K, Ca, Mg), nor the exchangeable acid cations Fe and Mn (data not shown), either for the first or the second year post-harvest.

Impacts of Harvesting on Fuel Load and Fire Risk Reduction
The pre-and post-harvest fire parameters for each of the sites are shown in Table 7. Before harvesting, site FA showed higher values for all the analyzed variables, indicating increased fire risk compared with NV or FB. This is due to the high average height of the broom plants in this area, as well as the presence of an abundant herbaceous load, resulting in greater fire propagation speed, heat per unit area, fire line intensity, and flame length than at the other sites ( Table 7). The effectiveness of shrub harvesting in reducing the risk of fire was greatest at site NV. Thus, considering an intermediate scenario with a wind speed of 15 km h −1 , the fire propagation speed was reduced by 84% compared with the pre-harvest scenario, 69% in terms of heat per unit area, 95% in terms of fire line intensity, and 74% in terms of flame length (Table 7). These reductions were less dramatic but also noteworthy for site FA (75%, 42%, 85%, and 59%, respectively, for the abovementioned parameters) and somewhat lower for site FB (57%, 40%, 74%, and 46% reductions, respectively).  Fire propagation velocity (S r ; m min −1 ), heat per unit area (H a ; kJ m −2 ), fire line intensity (F i ; kW m −1 ), and flame length (F l ; m). The percentage reductions in each variable with respect to the pre-harvest scenario are specified in parentheses.

Discussion
Mediterranean shrubs are frequently subjected to perturbations and disturbances associated with human activities. As a result, most of them are capable of regenerating easily and are highly resilient to these disturbances in the short term [13,25]. In particular, the rapid response of Erica australis and Calluna vulgaris following disturbance has previously been reported in the Iberian Peninsula [25], these species being capable of recovering their pre-disturbance spatial occupancy and cover values in a short space of time.
Our study reveals that all the locations studied presented significant decreases in general coverage one year after harvesting. However, two years after harvesting, there were no significant differences compared to the pre-harvest scenario at two of the three studied sites. The analysis of coverage at the level of the most important woody species revealed, however, that recovery in the second year after harvesting only occurred at sites with an abundant presence of ericaceous species. Ericaceous species were found to be those with the greatest sprouting capacity, as reported in other studies within the Iberian Peninsula [25]. The average height for all of the main shrubs species was, nevertheless, lower than in the pre-harvest scenario at all sites; therefore, the amount of biomass per hectare did not reach pre-harvest values.
The results obtained in this study reveal biodiversity index recovery rates of between 30 and 70% a year after harvesting, depending on the locality. These values are in line with those obtained in other studies which evaluated recovery rates of shrub formations within the Iberian Peninsula under very similar climatic conditions [25][26][27]. However, the values reached are far higher than those reported in other studies conducted in regions characterized by a colder climate, e.g., the French pre-alps [20], Canada, or the USA [28,29]. The latter studies reported that only between 15% and 20% of the species were able to recover after a period of one year following environmental disturbance. Other studies carried out under Mediterranean conditions involving shrub clearance in dehesa ecosystems in the Iberian Peninsula concluded that significant losses occurred in terms of species richness during the first years after clearance [9]. In addition, studies performed in very diverse shrub communities have reported that a minimum period of 3-5 years on average is required for real and complete regeneration of most shrub species after harvesting [25], and that in terms of productivity, a period of approximately 10 years is necessary to recover most of the shrub biomass extracted [19]. These aspects should be assessed at our study sites through future mid-term monitoring.
It has also been observed that the degree of dominance of some of the species increased after clearing at all the studied sites, especially one year after harvesting. In this regard, shrub species which were dominant prior to clearing have become dominated by herbaceous vegetation one year after the harvesting operation. This response is also in agreement with the general pattern observed in most other studies [25][26][27].
Mechanical shrub harvesting has also led to significant changes in the physical and chemical properties of soils at all the studied sites. Significant soil acidification has been observed in the three soil horizons of one location (NV), as well as in the deeper horizon of the other two locations (FA and FB) one year after harvesting. Some authors have previously reported the opposite effect (pH increase) [30,31], attributed to various factors such as humus transformation phenomena, release of cations by the decomposition of organic matter, or the use of hydrogen ions during the mineralization processes, among others. The machinery used for harvesting in our study left a large amount of wood on the soil (in some cases, losses over 60%), from which we hypothesized an increase in acid cations, despite the short period of time elapsed. This fact was verified though a cation exchange capacity analysis, which showed significant increases in aluminium and calcium depending on the location. Furthermore, the total litter thickness and the weight of the litter layer one year after harvesting had significantly increased in those locations and horizons where soil acidification was observed.
Harvesting also resulted in a significant increase in bulk density in some of the studied sites, both in the superficial horizon (of NV) and in the deep horizon (of FA). This effect was also observed in other studies carried out in the Mediterranean region which focused on tree species following forest harvesting operations [17] and was attributed to the process of soil compaction [32]. The increase in bulk density is also reflected by reduced soil porosity, restricting air movement and water permeability, according to [33]. In contrast, the opposite effect, i.e., reduction in bulk density, was observed at the FA site throughout the subsuperficial and deep horizons during the second year post-harvest. This effect may be associated with the increasing content of organic material at FA in the same soil horizons in which the bulk density was significantly lower. This negative correlation between organic material and bulk density seems to improve soil porosity [34,35].
We observed a significant increase in C and N after harvesting across the three studied soil horizons at two of the sites (FA and FB). This may be related to the fact that the degree of insolation reaching the soil increases considerably when the shrub mass is opened up by the clearing operation, so that the temperature of the soil tends to rise. Soil temperature is involved in many processes related to the rates of mineralization and decomposition of certain elements, including carbon and nitrogen [36], so the degree of mobilization of these elements could also have changed. The NV site, where we found no significant change in C and N contents, is located at a higher altitude than the others and temperatures are lower, so it is possible that all the mobilization and decomposition processes varied in comparison to the other two locations. In any case, these aspects should be assessed in future through more in-depth analyses.
Furthermore, the increase in C and N observed in deep horizons may be related with soil compaction leading to a migration of clay at depth. Clay minerals are involved in most of the organo-mineral associations of the soil; therefore, C and N contents increase with the clay content [37,38]. This phenomenon was verified in our study though a preliminary texture analysis of each horizon and has been previously been described in other studies [39].
A positive effect of harvesting in reducing the biomass load and, in turn, parameters linked to fire risk was found. Some studies of shrubland in the northwest of the Iberian Peninsula [40] have reported an average reduction of 50% in fire propagation speed 3 years after mechanical harvesting using the WFDS simulation model. Our results point to even higher percentage reductions in this variable at all the locations, reaching values of around 60% at FB and 83% at NV, although in this case, the results correspond to the situation 2 years after harvesting. As regards fire behavior (rate of spread, flame height, and fireline intensity), significant differences were found 2 years after clearing in a representative heathland in northwest Spain using a wind tunnel simulation experiment [13]. They found reductions of around 60% in the rate of spread and fireline intensity and 35% in flame height-findings which agree with those for FB in our study, but are lower than the values reached at NV or FA. This is probably due to the greater height of the shrubs at these sites prior to clearing.
In any case, it is also necessary to continue monitoring all these parameters over time in order to determine how long it takes for the resulting shrub formation to reach the pre-harvest fire risk level. The scarce studies undertaken in shrubland within the Mediterranean area point to an average timespan of approximately 4-5 years for a significant return to pre-harvest fire risk, depending on the type of shrub formation [13,41]. Hence, further research is needed to assess the temporal effects of this treatment on fire behavior in the mid to long term at the study sites.

Conclusions
Mechanical harvesting in shrubland growing under a temperate Mediterranean climate leads to positive and negative environmental impacts in the short term. These impacts should be estimated through a monitoring schedule for evaluating all the effects as a whole, i.e., considering the biodiversity and regeneration of species, soil properties, and fire risk at the same time to provide an overall picture for decision-making. A significant decrease in the total species coverage was found at all the locations one year after harvesting. However, the species coverage two years after harvesting was similar to the pre-harvest scenario at some sites due to the presence of ericaceous species. These ericaceous species displayed a high sprouting capacity, although their pre-harvest height was not recovered during the study period, or, therefore, their initial biomass. As regards biodiversity indices, the recovery rates after harvesting were very high, with values ranging from 30 to 70% depending on the location. Significant changes in the physical and chemical properties of soils were also observed. In this regard, slight acidification or soil compaction was found at some locations, as well as increments in carbon and nitrogen content after harvesting at others. Furthermore, two years after clearing, the fuel load and, therefore, the fire risk were found to be lower than in the pre-harvest situation at all the locations. However, mid-term impact assessments must be carried out to correctly evaluate possible changes in the vegetation, to estimate the rotation period of the shrub species, and to determine whether these changes are maintained over time or whether the shrubland returns to the pre-harvest situation given the plasticity that this type of ecosystem can display.