sustainability Article Responses in Soil Carbon and Nitrogen Fractionation after Prescribed Burning in the Montseny Biosphere Reserve (NE Iberian Peninsula) Sangita Chowdhury 1, *, José Manjón-Cabeza 1 , Mercedes Ibáñez 1,2 , Christian Mestre 1 , Maria José Broncano 3 , María Rosa Mosquera-Losada 4 , Josefina Plaixats 3 and M.-Teresa Sebastià 1,2 1 Group GAMES, Department HBJ, School of Agrifood and Forestry Science and Engineering (ETSEA), University of Lleida (UdL), 25198 Lleida, Spain; jose.manjon-cabeza@hbj.udl.cat (J.M.-C.); mercedes.ibanez@ctfc.es (M.I.); runge.mestre.christian@gmail.com (C.M.); teresa.sebastia@udl.cat (M.-T.S.) 2 Laboratory of Functional Ecology and Global Change (ECOFUN), Forest Science and Technology Centre of Catalonia (CTFC), 25280 Solsona, Spain 3 Grup de Recerca en Remugants, Departament de Ciència Animal i dels Aliments, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain; mariajose.broncano@uab.cat (M.J.B.); josefina.plaixats@uab.cat (J.P.) 4 Department of Crop Production and Engineering Projects, High Polytechnic School, University of Santiago de Compostela, 15782 Lugo, Spain; mrosa.mosquera.losada@usc.es * Correspondence: sangitafes@gmail.com Abstract: Prescribed fire is one of the most widely-used management tools to recover encroached rangelands. Fire has been reported to cause changes in the soil physical and chemical properties.  However, the legacy effects of former plant species on soil responses to fire remains unknown.  The legacy effect of the former extant plant species on soil carbon (C) and nitrogen (N) fraction- Citation: Chowdhury, S.; ation distribution after prescribed burning in topsoil (0–5 cm and 5–10 cm) was investigated in Manjón-Cabeza, J.; Ibáñez, M.; Mediterranean shrublands in Montseny. We sampled soils under five vegetation patch types: Cytisus Mestre, C.; Broncano, M.J.; Mosquera-Losada, M.R.; Plaixats, J.; scoparius L., Calluna vulgaris L., Erica arborea L., Pteridium aquilinum L., and Cladonia biocrusts, pre- Sebastià, M.-T. Responses in Soil and post-burning. Multivariate analysis on soil C and N fractions showed that soils under the legume Carbon and Nitrogen Fractionation Cytisus and the biocrust were the most differentiated. Vegetation patch types tended to respond after Prescribed Burning in the differently to burning, soils under Cytisus, Cladonia and Calluna showing the strongest response. Montseny Biosphere Reserve Total C and N, and C and N in sand decreased after burning in the 0–5 cm soil layer. Conversely, (NE Iberian Peninsula). Sustainability C in silt, as well as N in clay and silt, increased with soil depth after burning. This study will be 2022, 14, 4232. https://doi.org/ helpful for understanding ecological legacy effects and their possible consequences when planning 10.3390/su14074232 prescribed burning. Academic Editor: Ivo Machar Keywords: prescribed burning; soil particle size fractions; plant species-fire interactions; Cytisus Received: 18 November 2021 scoparius; Calluna vulgaris; biocrusts Accepted: 29 March 2022 Published: 2 April 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in 1. Introduction published maps and institutional affil- Fire is an important driver of environmental changes in an ecosystem, responsible for iations. altering nutrient pools by changing the physical, chemical and biological properties of soils and nutrient cycling [1,2]. Prescribed burning is considered as the deliberate application of fire characterized by lower temperature, intensity and severity compared to wildfires [3,4]. Copyright: © 2022 by the authors. This practice is a widely-used to reduce the risk of wildfire and recover shrublands [5]. Licensee MDPI, Basel, Switzerland. Soils are considered as a major reservoir of C in terrestrial ecosystems [6]. Soil includes This article is an open access article different carbon (C) and nitrogen (N) pools where soil labile C and N pools are character- distributed under the terms and ized by their small size and fast turnover rates, with large size and slow turnover rates, conditions of the Creative Commons recalcitrant organic C and N pools are physically or chemically protected [7]. Soil organic Attribution (CC BY) license (https:// carbon (SOC) is responsible for soil fertility, and it is used as an indicator of soil health. De- creativecommons.org/licenses/by/ creased soil C/N ratio led to a decline in SOC and an increase in soil N [8,9]. The influence 4.0/). of prescribed fires on soil C and N depends on fire frequency and intensity, in addition to Sustainability 2022, 14, 4232. https://doi.org/10.3390/su14074232 https://www.mdpi.com/journal/sustainability Sustainability 2022, 14, 4232 2 of 14 climatic factors [3,10]. Furthermore, prescribed burning is responsible to affect only the upper centimetres of the soil [4]. Prescribed fires may also affect microbial community composition and function that lead to altered C and N cycling in an ecosystem [10,11]. Some studies have described increased upper soil C and N content after low intensity fires because of the incorporation of unburned or partially unburned slash fragments into the soil [12–14], while other studies report no change [15]. Soil organic matter decreased with short-term low intensity fire in topsoil in Mediterranean Shrubland [16,17]. [18] observed decreases in soil nutrient concentrations with greater fire intensity on topsoil in a mixed Chaparral in USA. Soil carbon and nitrogen were lost from the first 5 cm with high intensity and long-term burning in a shrubland of Eastern Spain [19]. Soil particle size fractions are helpful in maintaining the stability of ecosystems [20]. Soil organic matter within the sand fraction is allocated to the active (labile) pool, and that in the silt and clay fractions, to the passive (recalcitrant) pool. The labile pool is easily affected by fluctuation in environmental conditions as well as decomposed rapidly and becomes oxidized easily with changes in land use [21]. The passive/non-labile pool is more stable and recalcitrant, and therefore this fraction is decomposed slowly by microbial activity [22]. Furthermore, plants determine the quantity and the quality of residues, soil organic matter, as well as soil structure [23]. Thus, soil functions are also affected by plant func- tional diversity [24]. Plant functional types (PFTs) have proved to be a useful tool for predicting soil processes including C and N cycles [25]. Legumes have the potential to modify soil nutrient availability as legumes have the capacity to fix symbiotic N [26]. In addition, biological soil crusts (BSCs), which are assemblages of lichens, fungi, cyanobac- teria, and mosses that colonize the soil surface, play a key role in the N cycle, because N-fixing lichens and free-living heterotrophic bacteria forming part of BSCs are able to fix substantial amounts of atmospheric N [27]. Many studies have focused on the dynamics of topsoil C and N stocks after prescribed burning. There are several studies conducted on fire effects on seed germination, natural regeneration, and changes in vegetation composition. However, very few studies have focused on the effects of fire on soil C and N fractionation, that is, C and N contents in clay, sand, and silt fractions [16,28]. Furthermore, to our knowledge, this is the first study that examines the legacy effects of the prior plant species on soil carbon and nitrogen particle size fractions after prescribed burning. Ecological legacy is a concept focused mainly on community or ecosystem-level phenomena, related to memories of the ecosystem to past events [29]. To understand present biodiversity patterns and predict the future ecological impacts of ongoing human practices on ecosystem services and functions, it is important to understand the carry-overs, or legacies, of the past events [30]. Therefore, the aim of this study was to assess the impact of prescribed burning on the C and N contents in the different fractions in topsoil (0−5 cm and 5−10 cm) in Pla de la Llacuna, Montseny, particularly to examine the legacy effect of the former extant plant species on soil carbon and nitrogen fractions after prescribed burning. We hypothesise that due to low to moderate prescribed fire, (1) the C and N contents in different soil fractions will be increased because of accumulation of burnt material; (2) patches dominated by different plant functional types will show variability in the soil C and N contents, with soil under legume containing higher C and N contents compared with other species; and (3) burning effects on soil C and N distribution will be modified by prior plant species, that is, there will be ecological legacy effects of former extant species. 2. Materials and Methods 2.1. Study Area The study was conducted at the Pla de la Llacuna (longitude 2◦ 180 to 2◦ 220 east, latitude 41◦ 440 to 41◦ 470 north), located in the Pla de la Calma, an elevated plateau that ranges between 1000 and 1350 m a. s. l. This occupies an irregular area of about 974 hectares in the Montseny Natural Park, in the Northern Catalan Pre-Littoral Range, north-eastern Iberian 2. Materials and Methods 2.1. Study Area The study was conducted at the Pla de la Llacuna (longitude 2°18′ to 2°22′ east, latitude 41°44′ to 41°47′ north), located in the Pla de la Calma, an elevated plateau that ranges between 1000 and 1350 m a. s. l. This occupies an irregular area of about 974 Sustainability 2022, 14, 4232 3 of 14 hectares in the Montseny Natural Park, in the Northern Catalan Pre-Littoral Range, north- eastern Iberian Peninsula (Figure 1). This park comprises Mediterranean and Central European landscapes including different biomes, having local influence of metropolitan Peninsula (Figure conurbations 1). This park comprises Mediterranean and Central European landscapes nearby. including different biomes, having local influence of metropolitan conurbations nearby. Figure 1. (A) Study area located in the NE of the Iberian Peninsula, (B,C) in the Natural Park of Figure 1. (A) Study area located in the NE of the Iberian Peninsula, (B,C) in the Natural Park of Montseny(Catalonia). Montseny (Catalonia).(D) (D)RGB RGB orthomosaics orthomosaics representing representingthe thepre-burned pre-burnedand andburned burnedarea of of area thethe Pla de lade Pla Llacuna experimental la Llacuna area (courtesy experimental of the area OPEN2PRESERVE (courtesy of the project, https://open2preserve. OPEN2PRESERVE project, eu/ accessed on 17 November https://open2preserve.eu/ 2021). accessed on 17 November 2021). The plateau The plateau isischaracterized characterizedby by humid Mediterranean humid climate; Mediterranean mean annual climate; mean precipita- annual tion approximately ranges between 700 and 1000 mm, where precipitation approximately ranges between 700 and 1000 mm, where snow accounts for around snow accounts for 10% of the annual ◦ around 10% of precipitation. The mean annual the annual precipitation. temperature The mean is 11 C according annual temperature is 11 °C to the Taga- according to the Tagamanent meteorological station (https://www.meteo.cat accessed on2021). manent meteorological station (https://www.meteo.cat accessed on 17 November 17 Bedrock is a metamorphic schist where the major minerals are quartz albite, muscovite November 2021). Bedrock is a metamorphic schist where the major minerals are quartz and chlorite. Soils are acidic, with a pH of 4.5 to 5.5, and characterised by a sandy-loam albite, muscovite and chlorite. Soils are acidic, with a pH of 4.5 to 5.5, and characterised texture [31]. The topographic location and humidity circulation lead to a general situa- by a sandy-loam texture [31]. The topographic location and humidity circulation lead to a tion that favors the presence of Atlantic vegetation, where the shrubs Erica scorparia L., general situation that favors the presence of Atlantic vegetation, where the shrubs Erica Erica arborea L. and Calluna vulgaris L. are widespread. The area is covered principally scorparia L., Erica arborea L. and Calluna vulgaris L. are widespread. The area is covered by shrubs and grass [32], including biological crusts dominated by various species of the principally by shrubs and grass [32], including biological crusts dominated by various lichen Cladonia sp pl., usually accompanied by some grasses. For centuries, the hills of the species of the lichen Cladonia sp pl., usually accompanied by some grasses. For centuries, Montseny mountains have been used as pastures. Montseny was declared as a Natural Park the hills of the Montseny mountains have been used as pastures. Montseny was declared in 1977; since then many traditional practices have declined, and some, such as shepherds’ as a Natural Park in 1977; since then many traditional practices have declined, and some, burning, have completely ceased. Many small flocks composed mainly of sheep and goats such as shepherds’ grazed in Montseny burning, have in former completely times. ceased. The lack Many of direct smallhas grazing flocks composed caused a changemainly in the ofvegetation sheep andstructure goats grazed in Montseny in former times. The lack and open grasslands are now covered by shrubs [33]. of direct grazing has caused a change in the vegetation structure and open grasslands are now covered by shrubs [33]. 2.2. Field Sampling and Laboratory Determinations The area selected for the prescribed burning experiment had a surface of 1.7 ha (Figure 1). A smaller area of 30 × 90 m within this surface was selected for the soil sampling (Figure 1). The burning event was conducted in 28 February 2019. Two samplings were performed in the study area, one before the burning event (pre-burning) which was carried out in 30 January 2019; and a re-sampling afterwards (post-burning) in 5 March 2019. The fire intensity was low to moderate in our study area. There was a small rain after burning in February 21 (0.1 mm). Sustainability 2022, 14, 4232 4 of 14 Five vegetation patch types were initially identified as the most abundant vegetation patches in the shrubland. Each type was dominated by different species: the vascular plants Calluna vulgaris L., Erica arborea L., Cytisus scoparius L. and Pteridium aquilinum L.; and the biocrust dominated by lichens mostly of the genus Cladonia. These vegetation patches belong to different plant functional types, where the shrubs C. vulgaris and E. arborea are Ericaceae and C. scoparius is a legume; P. aquilinum is a fern; and the Cladonia patch is a lichen- dominated biocrust, often mixed with grasses. Henceforth, we will use the terminology: Calluna (CV) for Calluna vulgaris patches; Erica (EA) for Erica arborea patches; Cytisus (CS) for Cytisus scoparius patches; Pteridium (PA) for Pteridium aquilinum; and Cladonia (CSP) for the biocrusts with Cladonia. The sampling was conducted using a stratified directed sampling. Soil samples were extracted in two different soil layers (0−5 cm and 5−10 cm) using a 4 × 4 cm2 coring probe. The treatments used for stratification (burning, vegetation patch and soil layer) were replicated six times across the sampling area, resulting in 120 samples. Sampling points were georeferenced throughout the sampling area with a highly precise GNSS Leica Zeno 20 (Leica Geosystems AG, Heerbrugg, Switzerland) with a differential correction Real Time Kinematic (RTK) broadcast system that was connected to the RTKAT service. In addition, we measured the slope at each soil sampling point using the 5 m DEM (Dig- ital Elevation Model) from the IGN (www.ign.es accessed on 17 November 2021). Slope was calculated in QGIS based on the DEM provided by the Institut Cartogràfic Geogràfic Català (ICGC, 2021), and a slope value provided for each sampling point. Slope was included in the modelling as a way of controlling the possible effects of spatial microtopographical heterogeneity. Afterwards soils were transported to the laboratory and oven dried at 60 ◦ C until constant weight. Soil samples were physically fractionated using the method developed by Six et al. [34]. The three soil fraction samples of known moisture content were analysed by a LECO C.N.H.S. Elemental Analyzer for the percentage of C and N. Total percentage of soil C and N was also measured in the two different soil depth layers. This results in eight variables including: total C, C in clay, C in silt, C in sand, as well as total N, N in clay, N in silt, N in sand. The C/N ratio was also calculated for each of the fractions described above. 2.3. Data Analysis Multivariate indirect ordination analysis was applied to the ensemble of the soil variables analysed in this study. In particular, we applied Principal Component Analysis (PCA). Multivariate analysis was conducted with CANOCO 5.1 [35]. In addition to multivariate analysis, we performed univariate statistical regression on each study variable, using linear mixed effect models, with the identity of the sampling point as random factor, and with burning, vegetation patch type, soil layers and slope as fixed factors. We included the interactions among those factors, except for slope, which was used as a covariate. Then the best model was selected according to AIC (Akaike Information Criterion) and stepwise method (both forward and backward). We did the normality test of all variables by Shapiro-Wilk Normality Test. All the variables followed the normal distribution except C/N ratio in the clay fraction. In this case, we used a generalized linear mixed model with inverse gaussian distribution. Then we performed Tukey’s post hoc analysis at the p < 0.05 significance level. All statistical analyses were conducted in the R version 4.0.2 [36]. 3. Results In the study area, silt was the fraction including the highest soil C and N contents (Tables 1 and 2). We found that most of the study variables increased with slope as well as soil depth (Tables 1–3). Sustainability 2022, 14, 4232 5 of 14 Table 1. Final models for all C variables in the study. p-values from the mixed model regressions on C variables, for the following tested explanatory variables: slope, burning, species, and soil depth. Explanatory Carbon Variables Variables Total Soil Carbon in Carbon in Carbon in Silt Carbon Clay Sand Slope 0.001 ** 0.001 ** 0.002 ** <0.001 *** Burning 0.004 ** 0.175 0.191 0.085 Species 0.095 0.061 0.076 0.024 * Soil depth <0.001 *** <0.001 *** <0.001 *** <0.001 *** Burning * Soil depth 0.021 * - 0.039 * 0.034 * Species * Soil depth - - - 0.076 p < 0.1; p < 0.05 *; p < 0.01 **; p < 0.001 ***. 3.1. Overall Soil C and N Distribution Considering together all the C and N variables, including C and N content in the total soil and in the three soil fractions, the Principal Component Analysis (PCA) explained 65% of the total variability. PCA axis 1 explained 88% of the explained variability and mainly separated samples according to the plant species under which the soil was originally extracted (Figure 2A). Soils under the legume shrub Cytisus were distributed along the most positive part on PCA1, followed by those under Pteridium (Figure 2A); while soils under Cladonia distributed mostly on the negative side (Figure 2A). Soils under the two Ericaceae had low responses to this axis, suggesting intermediate soil C and N trends (Figure 2A). In addition, PCA axis 2 added 6% to the explained variability. This axis mostly explained the variability due to burning (Figure 2A). Soils under Cytisus and Cladonia followed by Calluna were the ones showing a higher response to burning compared to the soil under the other species (Figure 2A). Finally, PCA axis 3 explained 3% of the total explained variability. Axes PCA2 and PCA3 suggest differences among vegetation patches in overall soil C and N responses to burning, soils under Cytisus, Calluna and Cladonia being the ones most differentiated between unburned and burned conditions (Figure 2A,B). 3.2. Total and Fractional Soil C Both total C (Figure 3A; Table 1) and C in the sand fraction in the 0–5 cm soil layer (Figure 3D; Table 1) decreased after burning; conversely, C in the silt fraction (Figure 3C; Table 1) increased with prescribed burning in the 5–10 cm layer, according to the regression model. However, this was not captured by the Tukey test, suggesting that the increase of C in the silt fraction with burning is relatively weak. The highest C content was found in the soils under Cytisus patches, while the lowest C content was found under Cladonia patches in the sand fraction (Figure 4D; Table 1). No significant interactions were found between plant species and burning in the total carbon, neither in the three soil C fractions (Table 1). However, we could detect tendencies in the particular response of some plant species to burning (Figure S1 in Supplemen- tary Materials). After burning, total soil C in the 0–5 cm soil layer showed a decreasing trend in all species, which was more pronounced in Cladonia and Erica patches compared to the other patch types (Figure S1A in Supplementary Materials). C in sand also de- creased in all the species, but especially in the Cytisus and Calluna patches in the 0–5 cm soil layer (Figure S1D in Supplementary Materials). In contrast, soil C increased after burning more remarkably in Cytisus and Calluna in the silt fraction in the 5–10 cm layer (Figure S1G in Supplementary Materials). Sustainability 2022, 13, x FOR PEER REVIEW 6 of 15 Sustainability 2022, 14, 4232 6 of 14 Samples Figure2.2.Samples Figure distribution distribution alongalong thefirst the two two first axes axesand (PCA1 (PCA1 PCA2and (A);PCA2 (A); and the firstand andthe first and third third axes,PCA1 axes, PCA1andand PCA3 PCA3 (B)(B) of Principal of Principal Component Component Analysis Analysis (PCA)(PCA) performed performed on the soil on the overall overall C soil C and and N variables N variables (total (total andand fractional)in fractional) in 0–5 0−5 and and 5−10 5–10soil soillayers. Samples layers. areare Samples clumped clumpedinto into groupsgroups accord- according to vegetation patch type and prescribed burning treatment. The mean value ± 1 standard ing to vegetation patch type and prescribed burning treatment. The mean value ± Sustainability 2022, 13, x FOR PEER REVIEW 1 standard 7 of 15 error error of each axis variable is represented for each species and treatment by different symbols and of each whiskers.axis variable is represented for each species and treatment by different symbols and whiskers. 3.2. Total and Fractional Soil C Both total C (Figure 3A; Table 1) and C in the sand fraction in the 0−5 cm soil layer (Figure 3D; Table 1) decreased after burning; conversely, C in the silt fraction (Figure 3C; Table 1) increased with prescribed burning in the 5−10 cm layer, according to the regression model. However, this was not captured by the Tukey test, suggesting that the increase of C in the silt fraction with burning is relatively weak. Figure 3. Mean ± 1 SE total and fractional (top to bottom) soil C (A–D), N (E–H) and C/N ratio (I–L) Figure 3. Mean ± 1 SE total and fractional (top to bottom) soil C (A–D), N (E–H) and C/N ratio (I–L) distribution distribution inin thethe 0−50–5cm cm andand 5−105–10 cmlayers. cm soil soil layers. Different Different letterssignificant letters indicate indicate significant differences differences (p << 0.05) (p 0.05)among among treatments treatments according according to multiple to multiple Tukeycomparison Tukey mean mean comparison tests (soil tests (soil layers and layers and burning). burning).Sample Sample sizesize n =n 120. = 120. The highest C content was found in the soils under Cytisus patches, while the lowest C content was found under Cladonia patches in the sand fraction (Figure 4D; Table 1). No significant interactions were found between plant species and burning in the total carbon, neither in the three soil C fractions (Table 1). However, we could detect tendencies in the Sustainability Sustainability2022, 2022,13, 14,x4232 FOR PEER REVIEW 7 of814of 15 Figure 4. Mean ± 1 SE total and fractional (top to bottom) soil C (A–D) and N (E–H) distribution per vegetation patch type; Cytisus scoparius (CS); Cladonia biocrust (CSP); Calluna vulgaris (CV); Erica arborea (EA); Figure 4. Pteridium Mean ± 1 SE aquilinum total and (PA). Different fractional letters soil (top to bottom) indicate significant C (A–D) differences and N (E–H) (p < 0.05) distribution amongFigure 4. Meanpatch pertreatments vegetation ± 1 SE totalCytisus according type; and to fractional Tukey mean scoparius (top to bottom) comparison (CS); Cladonia soil C (CSP); tests. biocrust (A–D) Sampleand N (E–H) size Calluna distribution n = 120. vulgaris (CV); Erica per vegetation patch arborea (EA); type; aquilinum Pteridium Cytisus scoparius (CS); Cladonia (PA). Different biocrustsignificant letters indicate (CSP); Calluna vulgaris differences (p (CV); < 0.05)Erica arborea (EA); Pteridium aquilinum (PA). Different letters indicate significant differences (p < 0.05) among Soil Ctreatments content according in the 5−10to Tukey cm mean layercomparison was lesstests. thanSample thann =that half size 120. in the 0−5 cm layer among treatments according to Tukey mean comparison tests. Sample size n = 120. under Cytisus, Soil C Erica contentand Pteridium. in the 5–10 cm While soilless layer was C content than halfinthan the that 5−10incm thelayer was 0–5 cm half than layer that in the under 0−5 Soil C cm Cytisus, layer content Erica inunder and the Cladonia 5−10 cmWhile Pteridium. and layer soil Calluna was (Figure less than C content in half 5; Table than the 5–10 that cm 1). in the layer was0−5 halfcm layer than that in Cytisus, under the 0–5 cm layer Erica andunder Cladonia Pteridium. While Calluna andsoil (Figurein5;the C content Table 1).cm layer was half than 5−10 that in the 0−5 cm layer under Cladonia and Calluna (Figure 5; Table 1). Figure 5. Mean ± 1 SE soil C distribution in sand fraction in 0−5 cm and 5−10 cm soil layers per Figure 5. Meanpatch vegetation ± 1 SE soil type; C distribution Cytisus in Cladonia scoparius (CS); sand fraction biocrustin(CSP); 0−5 cm andvulgaris Calluna 5−10 cm soilErica (CV); layers per Figure vegetation 5. Mean ± 1 patchPteridium arborea (EA); SE soil type; Cytisus C aquilinum distribution scoparius in sand (CS);letters (PA). Different fraction Cladonia in biocrust indicate 0−5 cm (CSP); significant and 5−10 Calluna differences cm soil (p