Effects of firewood harvesting intensity on biodiversity and ecosystem services in shrublands of northern Patagonia

Forest management has historically focused on provisioning of goods (e.g. timber, biomass), but there is an increasing interest to manage forests also to maintain biodiversity and to provide other ecosystem services (ES). We evaluated the effects of firewood harvesting intensity on biodiversity and different ES in three contrasting shrubland sites in northern Patagonia (Argentina). At each site, four harvesting treatments, representing various levels of harvest intensity, were randomly assigned to eight permanent sample plots of 31.5 m × 45 m during 2013–2014. We found that the effects of increasing harvesting intensity on plant diversity changed from negative to positive (and from nonlinear to more linear responses) with increasing site productivity. Harvesting intensity showed contrasting effects on variables related to fire protection ecosystem service, since it reduced fuel amount (potentially reducing fire spread) but also reduced live fuel moisture content (potentially increasing flammability) at the three sites. Two variables related to soil formation and protection ES, leaf litter cover and aerial soil cover, decreased with harvesting intensity at the three sites. We conclude that shrubland management for firewood production may enhance biodiversity without compromising certain important ES. The intensity of harvesting should be determined according to site conditions and forecasted impacts on biodiversity, fire and soil formation and protection.


Background
Historically, forest management has focused on favoring the sustainable production of few items with commercial value (e.g. timber, pulp, biomass for energy) whereas today the concept of sustainability also implies the maintenance of biodiversity and multiple ecosystem services (ES). However, the response of biodiversity and multiple ES to management is often difficult to predict because ES may respond in different and often nonlinear ways (Steffan-Dewenter et al. 2007;Bennett et al. 2009;Carpenter et al. 2009). Although the topic of ES has been popular in ecology research, experimental tests of relationships between natural resource management and multiple ES are still necessary (de Groot et al. 2010).
Forest management activities (e.g. thinning, coppicing, harvesting) generally enhance provisioning ES (e.g. timber, biomass) but could have negative effects on supporting ES such as soil formation and protection (de Groot et al. 2010;Cimon-Morin et al. 2013;Biber et al. 2015). This trade-off generates externalities that are not incorporated when decision making is driven by financial analyses focused on products. On the other hand, there could be win-win scenarios, considering that some forest interventions may enhance understory diversity and species richness in afforestation's by releasing resources at soil level (Battles et al. 2001). In this way, multiple ES response to management might be complex.
The effects of forest intervention intensity on plant diversity can vary according to site conditions. The Dynamic Equilibrium Model (DEM) proposed by Huston (2014) predicts that site productivity will modulate the effects of harvesting intensity (seen as a disturbance gradient) on diversity (Huston 2014). Sites with greater productivity should tolerate more (or even benefit from) harvesting intensity.
Biomass harvest is an effective practice to reduce fuel load and fire spread probability (Regos et al. 2016) and consequently could improve the fire protection ES that mitigates or prevents fire damage to human health, safety, and livehoods (e.g. the usage of natural resources, tourism) (Haines-Young and Potschin 2018; Sil et al. 2019). This is especially important in a changing climate where temperatures are expected to rise (Halofsky et al. 2017). Nevertheless, forest intervention changes site environmental conditions. Forest openings increases temperatures, fluctuations of temperatures and wind speed, and reduces relative humidity (Trentini et al. 2017). These environmental changes at stand level may favor a reduction in live fuel moisture content (LFMC) of foliage on the remnant stand (Pollet and Omi 2002). LFMC is a key trait in determining communities' flammability being correlated with several flammability traits (Cornelissen et al. 2003). Therefore, partial harvesting could increase fire spread probability through changes in LFMC content in the remnant stand but also reduce it by changing the fuel amount of the stand. This has been already documented at a global scale (Oddi 2018) and may be expected to occur at stand level through harvesting intensity gradients.
Soil formation and protection are one of the largest ES contributing to total ecological service value of land (Liu et al. 2012) and are largely associated to litter production and soil cover (Cotrufo et al. 2015). Partial harvesting or thinning reduces stand density affecting litter production and soil cover (Harrington and Edwards 1999;Huebschmann et al. 1999;Roig et al. 2005;Jandl et al. 2007). Such interventions may reduce the organic layer of the forest floor, compromise soil protection from extreme changes in moisture and temperature, reduce water infiltration capacity, and promote soil erosion due to lower soil cover (Osman 2013). In addition, litter fall is a particularly important process in nutrient cycling of forest ecosystems, as it provides the main above-ground contribution of carbon and nutrients to the forest floor through the formation of humus (Gallardo et al. 1998;Chapin et al. 2012). Therefore, soil formation and protection are a key ES to incorporate in forest decision making.
Firewood management of north Patagonia shrublands could be economically sustainable (Goldenberg et al. 2018) but the effect of management on biodiversity and ES is yet unknown. In this sense, National Law 26.331 rules that extractive activities in native forests should follow a sustainable silvicultural plan that secures the long-term provision of forest related ES. In this context, local information on forest management is a key aspect for natural resources management (Rusch et al. 2017).
The objective of this study is to analyze the response of biodiversity and two ES (fire protection and soil formation and protection) to a harvesting intensity gradient under contrasting site conditions in northern Patagonian shrublands. The specific hypotheses are: plant diversity in sites with greater productivity would benefit from a higher harvesting intensity; increasing harvesting intensity reduces LFMC content in the remnant stand but also reduces the fuel amount of the stand and leaf litter cover and aerial soil cover are reduced by increasing harvesting intensity.

Study area and experimental sites
In the Andean region of north Patagonia, Argentina, there is a large area covered by pure and mixed Nothofagus antarctica shrublands. These forests types are dominated by resprouting woody species being the most diverse ecosystem of the region (Gyenge et al. 2009;Speziale et al. 2010). Usually, in absence of fire and cattle, they are displaced by longer-lived species and obligate seed-dispersers such as Austrocedrus chilensis or Nothofagus dombeyi trees forming pure forests (Kitzberger and Veblen 1999). Nonetheless, this dynamic usually is affected by fire due to fire-prone characteristics of these communities (i.e. fire adapted forests) (Keeley et al. 2011;Blackhall et al. 2012).
The north Patagonian Andean region is a Mediterranean type-climate region (Keeley et al. 2011), with an annual precipitation ranging from 920 to 1300 mm in the shrublands distribution that occurs mainly during autumn and winter season (Reque et al. 2007). Average annual temperatures range from 7°C to 9°C, with annual average maximum temperature of 15°C and minimum temperature of 1.5°C. Frosts occur about 120 days a year, with 0.5 days' hail, annual relative humidity 65%, and an annual dew temperature of 2°C (Reque et al. 2007;Gyenge et al. 2009). Presence of frost is longer in valley bottom, and cold air accumulation can limit tree growth (Tejera and Davel 2004).
Three sites with contrasting environmental conditions in the province of Rio Negro, Argentina were chosen to conduct the study (Fig. 1). The choice was based on site aspect, one of the main environmental factors driving shrubland physiognomy in this region. In south-aspect hillsides, soils are deeper, and have higher moisture retention than north-aspects hillsides, where soils are drier because they are exposed to the dominant northwestern winds and intense summer droughts (Davel and Ortega 2003). Consequently, productivity, in terms of cubic meters of firewood production (mean annual increment; Table 1), is determined by the position of sites in the landscape. Hence, high productivity site was located on a southern aspect, medium productivity site on a northern slop and low productivity site on a valley bottom where environmental conditions are limiting mainly due to cold air accumulation, high annual thermal amplitude and lower precipitations which favor the intense summer drought (Table 1). Overstory vegetation was dominated by mixed N. antarctica shrublands, Schinus patagonicus, Lomatia hirsuta, and Embothrium coccineum codominated the stands in the high and medium productivity sites, and pure N. antarctica shrublands at the low productivity site (Table 1); the only woody species in all three sites was N. antarctica. In the high and medium productivity sites, soils were dominated by the group of Hapludands, with dark color, sandy texture, poor structure and abundant presence of roots. The medium productivity site had shallower soils and presence of rocks. The groups of Udivitrands were dominant in the low productivity site with ocher color and poor abundance of roots. Soil chemical characteristics were similar between medium and low productivity sites but both poorer than the high productivity site (Carron et al. 2020). Elevation across sites ranges from 790 to 840 m.

Harvesting treatments
At each site eight 31.5 m × 45 m plots were selected (total plots = 24). Between 2013 and 2014 six plots were harvested in six strips of increasing width (1.5, 2.5 and 3.5 m, along the plots) with the two remaining plots serving as controls, resulting in approximately 0, 30%, 50% and 70% of shrub cover removal (i.e. harvesting intensity) respectively. All stems with more than 4 cm of diameter were harvested and classified as firewood leaving smaller branches and leaves in the intervention strips. The very few A. chilensis individuals that were present were not harvested. All firewood volume was estimated (Table 1).

Ecosystem services framework
We measured a suite of indicators (proxies) associated with plant diversity and key ES in each of the 24 plots between spring of 2015 and autumn of 2017. ES were classified as shown in Table 2. For this study, we limited our focus on: Biodiversity, fire protection and soil formation and protection.

Biodiversity: Shannon index, plant cover, richness
Plant diversity was determined by vegetation samples containing all plant life forms. In all plots each species cover was determined in quadrants of 1 m 2 . We used four quadrants per plot which were placed in the four cardinal points at 2.5 m from the plot center. Measurements were done in spring of 2015, autumn of 2016, spring of 2016 and autumn of 2017. In each measurement species were identified. We determined two indicators associated with vegetation: species richness (that is, number of species: S′) and plant cover (%). Thereafter, the Shannon index (diversity: H′) was calculated for plant diversity determination (Spellerberg and Fedor 2003): where P B is the relative abundance of each species.
Fire protection: live fuel moisture content, fuel amount We chose two indicators to evaluate the effect of shrubland management on fire protection ES: fuel amount and live fuel moisture content. Live fuel moisture content (LFMC) was measured during summer of 2016 and 2017 in all plots. Samples of live leaves and fine fuels were taken by crossing each plot perpendicular to harvesting strips direction. Mixed species samples (~100 g) were taken from each plot of the high and medium productivity sites and pure samples of N. antarctica (~50 g) on every plot of the three sites.
Leaves were harvested at 2-m height on north aspect (Bianchi and Defossé 2015). The compositions of the mixed samples (i.e. bags) were obtained considering the woody species composition of each plot (% cover area). The bags were stored inside fridges in field and then moisture content was obtained through gravimetric method (78 h at 62°C) in the laboratory (Cornelissen et al. 2003). This procedure was repeated every 2 weeks during the summer period (January-March), when fire risk is highest (Kitzberger et al. 1997). Fuel amount was determined using subplots in the center of each plot of 8 m of diameter. In each subplot we measured circumference at base height of all stems (approximately at 3 cm of ground level). Base diameter was then determined and the sum of the cross-sectional Table 2 Biodiversity, ES categories and indicators used in this work a a each plot represents the expected response of indicators to harvesting intensity. The hypothesis were mentioned in the introduction section with their respective references area of all trees was calculated as hectare surface (m 2 ·ha − 1 ).

Soil formation and protection: leaf litter cover, aerial soil cover
As an indicator of soil formation and protection, we used leaf litter cover and aerial soil cover. Aerial soil cover was calculated adding the remnant canopy cover (m 2 ·ha − 1 ) after harvesting to the understory cover measured with the vegetation samples as m 2 ·ha − 1 . Also, leaf litter cover was measured in 1-m 2 quadrants (four repetitions per plot) during autumn and spring in 2015-2016 and 2016-2017.

Data analysis
To analyze the effect of harvesting and site condition on each indicator, linear mixed-effects models were fitted using the lme function from nlme package in R 3.5 (R Core Team 2017) and glmer function and lme4 package when data follow non-normal distributions (Bates et al. 2014;Pinheiro et al. 2018). At the plot level, the models of plant cover, richness and Shannon index considered the fixed effect of harvesting intensity (quantitative predictor), squared harvesting intensity (to consider nonlinear responses to harvesting), site (categorical predictor), year (categorical predictor), season (categorical predictor) and interactions. Leaf litter cover model did not include the squared term as we expected a linear response (Table 2). Fuel amount and aerial soil cover considered the fixed effect of harvesting intensity (quantitative predictor), site (categorical predictor) and interactions. Finally, LFMC considered the fixed effect of harvesting intensity (quantitative predictor), date (categorical predictor including the four summer dates analyzed) and interactions. Plot was considered a random effect for avoiding pseudo-replications for plant cover, richness, Shannon index, leaf litter and LFMC. When necessary, variances were modeled using VarIdent function (for more details see Zuur et al. 2009;Pinheiro et al. 2018). Multi model inference was performed and AIC criteria selected the best models following a parsimonious criterion (Garibaldi et al. 2017) of all possible combinations using the dredge function and MuMin package (Bartoń 2009). The assumptions of homoscedasticity and normality were verified by visual evaluation of the residual scatter plots (e.g. residual vs. predicted values, q-q plots) and Shapiro-Wilk and Kolmogorov-Smirnov tests (shapiro.test and ks.test).

Biodiversity
The Shannon diversity index had a nonlinear response to harvesting intensity (Fig. 2). Nonetheless, in the selected model by AIC criteria (Table S1), the effect of harvesting intensity interacted with site. In the high productivity site, 70% harvesting enhanced diversity and reached the highest value of H′ for the site. The pattern was the same for both years and seasons (Fig. 2). In the medium productivity site, maximum diversity was achieved at intermediate harvesting intensity (50%) while in the lowest productivity site, maximum diversity was achieved in the plots with the lowest levels of harvesting intensity (0-30%). Shannon index was higher in the low productivity site when comparing the control plots (0% of harvesting intensity). Diversity was in general higher in the second year of measurements and the lowest levels were found for the first autumn. The magnitude of the harvesting effect in the low productivity site was less intense in spring (Fig. 2, Table S1).
Plant cover had nonlinear responses to harvesting intensity (Fig. S1) in high and medium productivity sites where intermediate harvesting enhanced the vegetation cover. The harvesting intensity effect was less clear in the low productivity site, where interactions with season were found. Plant cover usually was higher during the second year of sampling, except in the low productivity site where it was higher in autumn of the second year. When comparing among sites, the low productivity site had the highest levels of plant cover for the control plots (Fig. S1). Richness also had a nonlinear response to harvesting intensity although it strongly differed between sites (Fig. S2). In the high and medium productivity sites, intermediate harvesting enhanced species richness, which reached its highest levels at 50% of harvesting intensity at the high productivity site, and at 30%-50% in the medium productivity site. At the low productivity site, species richness was highest in the control plots and in low levels of harvesting intensity (0-30%). The negative effect of harvesting intensity was more marked in autumn than in spring. In this site, the only species in overstory is N. antarctica therefore the elevated levels of richness are mainly due to a rich understory composition. In the medium and high productivity sites overstory is more diverse as we found Lomatia hirsuta, Schinus patagonicus, Embothrium coccineum, Diostea juncea, and with less frequency Maytenus boaria, Discaria chacaye accompanying N. antarctica.

Fire protection ES
The general tendency was that LFMC decreased slightly as harvesting intensity increased (Fig. 3). N. antarctica had lower levels of moisture content than the mix vegetation in all the site conditions (Fig. 3). When comparing among sites, the lowest levels of moisture were found in the medium and low productivity sites. LFMC tends to decrease during the fire season independently of harvesting level, especially in mix vegetation (Table S2). Harvesting intensity reduced the fuel amount measured as the basal area of woody species (Fig. S3) and this effect was independent of site when fuel amount was measured as basal area (i.e. per ground area unit, m 2 ·ha − 1 ). Fuel Fig. 2 Harvesting intensity effect on Shannon index for three sites in two different years and seasons Goldenberg et al. Forest Ecosystems (2020) 7:47 amount was higher in the high productivity site than in the medium and low productivity sites (Fig. S3).

Soil formation and protection ES
Harvesting intensity had a strong negative linear effect on leaf litter cover (Fig. 4). In general, the high productivity site had higher levels of leaf litter and the low productivity site, the lowest. The higher levels of leaf litter were measured in the first year except in low productivity site. The lowest estimation of leaf litter cover was around 5% and was found in the low productivity site in plots treated with extreme harvesting intensity. The effect of harvesting intensity was lower in the second year (Fig. 4, Table S1). Aerial soil cover, composed by shrub canopy and herbs, decreased linearly with harvesting intensity at all sites (Fig. S4). The decrease in aerial soil cover was independent of the site effect (i.e. no interaction) when aerial soil cover was measured as basal area. Aerial soil cover  Goldenberg et al. Forest Ecosystems (2020) 7:47 was higher than hectare surface (> 10,000 m 2 ) because different strata were incorporated (shrub canopy and understory vegetation). Nonetheless, total soil cover (i.e. soil protection) was in general high if we consider the sum of the leaf litter cover (Fig. 4) to aerial soil cover (Fig. S4) in high and medium productivity sites.

Biodiversity
Our results suggest that shrubland management influences plant diversity in the northern Patagonia region. Based on our study which encompassed three different site conditions, we found that harvesting intensity affects plant species cover, richness and, consequently, diversity Shannon index, showing a strong interaction with site productivity. Because of the strong differences in site conditions, "site" was the most important predictor explaining plant diversity variability. We used Shannon index, that is one of the most applied diversity indicators in ecology and which resumes the interaction between plant cover and richness (Spellerberg and Fedor 2003). Following the DEM framework and the intermediate disturbance hypothesis low diversity would be expected for unmanaged and undisturbed sites, where woody species exclude understory species (Connell 1978;Battles et al. 2001), which is supported by our results in the medium productivity site. Similar results were found in other studies where high levels of woody trees diminished understory diversity due to competition for light in similar forest types (Quinteros et al. 2010) and other natural forests and plantations (Battles et al. 2001;Ares et al. 2010;Trentini et al. 2017). In this sense, Ishii et al. (2008) also found that linear thinning (similar to our strips harvesting) in Cryptomeria japonica plantations enhanced diversity, mainly due to changes in micro-climatic conditions. Ishii et al. (2008) showed that gap creation benefited vascular plants in thinned stands, although they did not explore an intervention intensity gradient (Ishii et al. 2008). In contrast, environmental conditions left after severe disturbances, such as those following heavy strip harvesting (low moisture, high deficit pressure vapor), could create stressful conditions suitable only for a very few species (Connell 1978). Therefore, diversity is expected to be maximized following medium intensity harvesting (50%), such as stated in the intermediate disturbance hypothesis (Connell 1978).
There are very few works that evaluated forest diversity through a disturbance gradient in South America temperate forests. One of these was carried out by Peri et al. (2016) in N. antarctica forests located in southern Patagonia, where this species forms tall forests due to site conditions and found that open forest favors species richness. Rusch et al. (2005) studying vegetation changes along a disturbance gradient in Chubut Province (south of Rio Negro, Fig. 1), have observed that in open stands (canopy cover~50%) understory cover is mainly dominated by adventitious grasses and herbs. Thus, forest interventions that partially opened the canopy resulted in an increase in plant richness and diversity supporting the intermediate disturbance hypothesis. In N. antarctica forests similar to those of our study (height < 8 m), Quinteros et al. (2010) when analyzing different stand structures also found that open stands (0-40% forest cover) enhanced richness and Shannon index of vegetation. Nonetheless, neither of these works had explored different site productivity conditions.
Our work demonstrates that the response of plant diversity to harvesting intensity interacts with site productivity, supporting the dynamic equilibrium model proposed by Huston (2014). Under this model, intermediate disturbance hypothesis is only expected to occur at intermediate site productivity (Huston 2014). Specifically, in our work, three patterns of response to increasing disturbance were found. In the high productivity site, extreme harvesting enhanced diversity probably due to intense disturbance that enhances light availability and counters a rapid occurrence of competitive exclusion by the dominant woody species (Huston 2014). Indeed, it is likely that at this site, high intensity harvesting does not create conditions too stressful for plants, due to the natural conditions of the site (high levels of soil moisture, attenuated extreme temperatures on summer). The intermediate productivity site showed the unimodal pattern proposed by Connell (1978) and other works in the region previously mentioned. In the low productivity site, extreme conditions leading to high plant mortality are achieved at low levels of disturbance and thus no/or very low levels of harvesting promoted diversity (Huston 2014). This is probably because water availability in summer is the most limiting resource in this site, so lower levels of soil water under high harvesting intensities decreases species richness (Stevens et al. 2006). In this site, spring weather conditions are moderate and thus the observed pattern during this season was more similar to the other sites. The most extreme conditions are reached in summer so autumn diversity shows how mortality of understory species is favored by higher levels of harvesting.

Fire protection ES
Fire protection ES is one of the main ES to be considered in decision making in northern Patagonia shrublands, since this fire-prone community is largely associated with urban centers which represents a danger to people and goods (de Torres Curth et al. 2012). The reduction of fuel load is one of the most widespread management practices to prevent fires in Mediterranean forests (Regos et al. 2016). However, the compromise between the effects of strips harvesting on fuel amount and moisture under continuous cover system has not been widely explored yet. Here we found that fuel amount, measured as basal area, decreased with harvesting intensity. Since fire is a contagious process (in a fire, each fuel particle is a source of ignition for the surrounding fuel particles) (Peterson 2002) the decrease in fuel amount could reduce the spread of fire. On the other hand, in all the sites we found that LFMC decreased linearly with harvesting intensity, in both mixed samples, which represent the specific composition of each site, as well as in those of N. antarctica, the control species found in all the sites. Therefore, harvesting could increase the shrubland flammability regardless of the quality of the site and the type of vegetation. Considering these two flammability traits there is a general trade-off between them.
The LFMC decreased throughout the summer for the entire harvesting intensity gradient, as might be expected for Mediterranean-type climate ecosystems (Pellizzaro et al. 2007). In these ecosystems, rainfall occurs in autumn-winter and water stress increases as the summer period passes. Although many species have strategies to keep their moisture level relatively constant throughout the year, the LFMC of N. antarctica behaves as seasonal (Bianchi and Defossé 2015). Regardless of the context, N. antarctica showed a lower LFMC than the mixed vegetation, meaning that it is one of the most flammable species of these communities (Ghermandi et al. 2016). In this sense, the low productivity site, formed by pure N. antarctica and high shrub cover due to extreme stems density of small dimensions (Table 1), demonstrates that shrublands under these conditions are highly fire-prone in accordance with previous works (Blackhall et al. 2017;Tiribelli et al. 2018).
In a context of climate change, it is key to understand the magnitude of the effect that management practices, such as partial harvesting, have on fuel amount and moisture, two of the most important drivers of fire regimes (Cornelissen et al. 2003). Given the trade-off between fuel amount and moisture, our results support the hypothesis that there is an intermediate harvesting intensity that would be optimal from the point of view of fire control. However, the effect size in LFMC was much lower than in basal fuel amount (Table S2). This would suggest that the focus should be on diminishing fuel loads independently of LFMC. Nonetheless, important components of fire behavior such as dead fuel loads (e.g. dead fine branches) and herbaceous loads, that can alter fire behavior, were not considered. Future studies should quantify the loads (t·ha − 1 ) and moisture in all of these components to improve quantitative evaluation of fire behavior (e.g. rate of spread, flame length) that would allow a deeper understanding of the multiple effects of harvesting intensity on the fire protection ES. Also, to determine the optimal harvesting intensity to reduce fire risk, future studies are required for assessing the magnitude of the continuity and moisture of the live fuels which modify the fire hazard of these communities.

Soil formation and protection ES
Leaf litter is one of the main detritus cover type in this community (de Paz et al. 2013). It is important for protecting the soil from erosion, increase aeration and regulate temperatures (Sayer 2006). In north Patagonia shrublands, leaf litter is mainly coming from woody species (de Paz et al. 2013). Our results ( Fig. 4) showed that harvesting intensity decreases leaf litter cover in the short term. This could affect soil properties having negative consequences in the long term considering that leaf litter is one of the main carbon inputs on the forest floors and a key component on soil formation process (Vesterdal et al. 1995;Jandl et al. 2007) therefore the soil formation and protection ES could be compromised. Our results show that leaf litter accumulation decreases with harvesting intensity, but the magnitude is diminished with time. Being that almost all the woody species can resprout (Rusch et al. 2017), leaf litter cover would increase as the sprouts get bigger which is observed in our experiment in the differences between treatments that tend to disappear in the second year of measurements. Probably, as woody species regrowth and the detritus input increases, the leaf litter cover would be reestablished as well as the aerial soil cover in the medium term.
Differences among sites in leaf litter cover are also crucial. The high productivity site had the highest levels of litter cover while the low productivity site had the lowest. In the low productivity site, regrowth is not as vigorous as in the intermediate and high productive sites (Goldenberg et al. unpublished work). All things considered, the low productivity site would be the most fragile in terms of soil formation and protection ES to the high levels of harvesting intensity (i.e. 70%). Nonetheless, to improve quantitative evaluation of soil erosion risk, soil erosion (t·ha − 1 ) should be directly quantified based on slope, soil and vegetation cover, under different harvesting scenarios that would allow a deeper understanding of the effects of harvesting intensity, on soil erosion risk and therefore the soil formation and protection ES.

Conclusions
In this work we show that biomass management in northern Patagonia shrubland affects biodiversity and different ES. In the high and medium productivity sites, biomass harvesting in strips enhanced biodiversity and did not show strong short-term negative effects on two important ES. The strongest negative effects of harvesting were found in the shrubland in the worst site conditions (i.e. valley bottom with extreme conditions) where increasing harvesting intensity had negative effects on biodiversity and soil formation and protection ES. Therefore, the intensity of harvesting should be determined according to site conditions and forecasted impacts on biodiversity, fire and soil formation and protection.