Habitat fragmentation is linked to cascading effects on soil functioning and CO₂ emissions in Mediterranean holm-oak-forests

Study area

The study area was located near Quintanar de la Orden (39°30′–39°35′N, 02°47′–02°59′W; 870 a.s.l.), in Toledo, southeastern Spain. This area has a Mesomediterranean climate characterized by 434 mm of mean annual precipitation and 14 °C of mean annual temperature (Ninyerola, Pons & Roure, 2005), with a pronounced summer drought, usually lasting from July to September. The landscape, a former predominant holm oak Mediterranean forest, is currently dominated by cereal and legume croplands, with scattered grape crops that complete the mosaic. The original forests are now highly fragmented in a variety of patch sizes, covering only the 28% of their original extent (Díaz & Alonso, 2003). The dominant tree is the holm oak (Quercus ilex L. ssp. ballota (Desf.) Samp), while the understory is mainly composed by shrubs of Kermes oak (Quercus coccifera L.) and, to a lesser extent, by species of Genista, Asparagus, and Rhamnus (for a full description of the study area see: Díaz & Alonso, 2003; Flores-Rentería et al., 2015; Flores-Rentería et al., 2016).

Experimental design and sampling

A total of three large (>10 ha; with 121 stems per ha on average) and five small (<0.5 ha; with at least three trees) forest fragments within an area of 1,000 ha, separated by a minimum of 50 m (to avoid spatial dependence) to a maximum of 8 km, were studied. Prevalent soils were Cambisols (calcics) (WRB, 2007), with sandy loam texture (17–39–44%, sand-silt-clay, respectively).

There is strong evidence that the exposure of the edges of fragmented forest causes changes in the abiotic and biotic conditions in comparison with forest interiors (Fischer & Lindenmayer, 2007; Flores-Rentería et al., 2015; Flores-Rentería et al., 2016; Murcia, 1995), while small fragments consist effectively only of edge habitat (Young & Mitchell, 1994). For instance, in the study area, forest interiors show higher intraspecific competition (measured as the proportion of area within a radius of 20 m from focal trees covered by other canopies), in comparison with edges and small fragments (0.46  ± 0.04, 0.36  ± 0.03, and 0.27  ± 0.14, interior, edges and small fragments, respectively) (Morán-López et al., 2016). Thus, we defined the influence of the agricultural matrix on forest fragments by the factor “matrix influence” with three levels: (1) low influence, at the interior of large fragments (at least 30 m from the forest edge; coded as “forest interior”); (2) mid influence, at the edges of large fragments (coded as “forest edge”); and (3) high influence, in small fragments (coded as “small fragments”). For each of the three large fragments, we selected five holm oak trees in the forest interior and five trees at the forest edge, while for each of the five small fragments we selected three holm oak trees (15 trees per fragmentation level), resulting in a total of 45 selected trees. For each of the focal trees, two coverage-sampling points were established: one under canopy (half way to the stem) and the other in open areas (with visible clear sky above, at least 1.5 m to another canopy), resulting in a total of 90 soil samples (see Flores-Rentería et al., 2016).

Field measurements

The field campaign was conducted in spring 2013, during the rainy growing season in Mediterranean ecosystems, when temperature and moisture were not limiting factors neither for plant growth nor for soil microbial functioning. Daily maximum temperature during the sampling days was of 24  ± 1.5  °C, with 15  ± 1.5 °C in average, and no precipitation was recorded.

Soil respiration (R_s) was measured at each sampling point with a portable dynamic closed chamber (SRC-1; P Systems, Amesbury, MA, USA) connected to an infrared gas analyzer (EGM-4; P Systems, Amesbury, MA). The soil chamber was placed on 90 external PVC collars (5 cm depth × 10 cm diameter) temporarily inserted to a depth of 3 cm into the soil 24 h before R_s measurements to minimize the impact of its insertion (e.g., underestimate R_s from roots) (Heinemeyer et al., 2011). R_s measurements were carried out at maximum daily soil activity (13:30–16:30 h) (Matías, Castro & Zamora, 2011), during three consecutive days (each sample point was measured one time). Immediately after R_s was measured, soil temperature and soil moisture were recorded at 10 cm depth by using a wireless multilogger thermometer (OMEGA, Norwalk, CT, USA), and a time domain reflectometer (TDR 300, Spectrum technologies, Illinois, USA), respectively. Then, soil cores (2 cm in diameter) were taken from a depth of 0–15 cm. Soils were sieved (<2 mm) and stored at 4 °C for later analyses.

A tree influence index was calculated for each sampling point, to take into account the influence of tree size and proximity according to the following formula: ${\displaystyle \text{Tree influence index}\left(\mathrm{Tii}\right)=\frac{\mathrm{BS}}{\mathrm{DT}}.}$

Where BS stands for the basal area measured at 25 cm from the ground (D₂₅), and based on stems ≥3 cm of diameter; and DT, the distance from the trunk, which was established as half the radio of the canopy cover in the case of the under-canopy points, and at 1.5 m from the edge of the canopy in the case of the open-areas. Additionally, height and canopy projection were measured for each holm oak tree. Tree basal area was selected to calculate the index, given its recognized direct relationship with soil functioning (Barba et al., 2013; Søe & Buchmann, 2005).

Soil properties

Soil moisture was determined by weight losses of 20 g samples oven-dried at 105 °C for 48 h. Total carbon (C) and nitrogen (N) contents were measured on air-dried soil samples, using a C:N elemental analyzer (Flash EA 1112 Series; Thermo Fisher Scientific, Waltham, MA, USA). Soil organic matter (SOM) was assessed by loss on ignition at 400 °C for 4 h (Ball, 1964). Microbial biomass C content was determined by the chloroform fumigation - extraction method modified by Gregorich et al. (1990).

Soil microbial community richness

The structure of soil bacterial and fungal communities was assessed by the DNA community fingerprinting technique of denaturing gradient gel electrophoresis (DGGE). The universal primers 338F/518R were used for amplification of the bacterial 16S rRNA gene (Muyzer, De Waal & Uitterlinden, 1993). In the case of fungi, the internal transcribed spacer nrDNA region ITS-1 was PCR-amplified using the primer pair ITS1-F/ITS2 (Gardes & Bruns, 1993). A GC clamp was respectively added to the 5′ end of forward bacterial (338F) and fungal (ITS1-F) primers to stabilize the melting behavior of the DNA fragments (Muyzer, De Waal & Uitterlinden, 1993). DGGE was carried out in gradients of 10–50% for fungi (Anderson, Campbell & Prosser, 2003) and 30–60% for bacteria (Grossman et al., 2010), with the concentrations of 7 M urea and 40% formamide (v/v) for the 100% denaturant. Electrophoreses were run at 60 °C 75 V for 16 h. Each band of the DGGE profile was hereafter referred to as an operational taxonomic unit (OTU). Gel bands were analyzed by using internal reference bands and known reference markers loaded in lanes at either side of the gel. The number of bands in a particular sample was considered comparative proxies of richness (S) of fungal or bacterial OTUs, respectively (Cleary et al., 2012). Similar analysis of DGGE banding patterns have been previously used in other studies (Anderson, Ellingsen & McArdle, 2006; Cleary et al., 2012; Flores-Rentería et al., 2015; Flores-Rentería et al., 2016; Suzuki et al., 2012; Vaz-Moreira et al., 2013).

Soil enzymatic activity

To characterize soil heterotrophic metabolism, we determined the polysaccharide-specific hydrolytic enzymes β-glucosidase and chitinase, and the phosphorus-mineralizing acid phosphatase in fresh soil samples. Enzymatic assays were based on methylumbelliferone (MU) (fluorogenic substrate) release upon cleavage by enzymes (Mathieu et al., 2013): MU-β-d-glucopyranoside (MU-G) for β-glucosidase (EC 3.2.1.3), MU-N-acetyl- β-glucosaminide (MU-Q) for chitinase (EC 3.2.1.14), and MU-phosphate free acid (MU-P) for acid phosphatase (EC 3.1.3.2). Stocks of 5 mM enzyme solutions were prepared in methoxy-ethanol, and enzyme substrates were diluted in sterile ultra-pure water to a final concentration of 800 µM for MU-P and 500 µM for MU-G and MU-Q assays. Stocks and calibration solutions, as well as diluted substrates were kept at −20 °C in the dark. Fluorogenic assays were performed by mixing 200 µl of soil supernatant (soil incubated over night with Tris-acetate buffer 10 mM, pH 4.5 in a horizontal shaker at 25 °C and 100 rpm), and 50 µl of the corresponding substrate with a final volume of 250 µl per well, using black 96-well microplates. Controls with soil supernatant heated at 100 °C for 10 min were also conducted separately for each sample. All reactions were performed at room temperature, applying a stirring of 500 rpm, in the dark and at different incubation times depending on the enzymatic test: 15 min (MU-P), 40 min (MU-G), and 60 min (MU-Q). After incubation, microplates were spin (3,000 rpm for 3 min), and 100 µl of the reaction mix was added to 100 µl of stopping buffer (Tris 1 M, pH 10–11). Measurements were carried out with a Victor3 microplate reader (Perkin-Elmer Life Sciences, Waltham, MA, USA), at 355/460 nm excitation-emission wavelengths. Experimental calibrations of known MU concentrations were performed, allowing estimating each enzymatic activity by extrapolating well fluorescence signals on the respective calibration regression lines. Blanks with buffer and fluorogenic substrates related to auto-fluorescence, and controls were subtracted from all measures. All enzymatic activities are expressed in pmol min⁻¹ mg⁻¹ of dry soil.

Data analysis

Prior to statistical analyses, all variables were tested for normality, and log transformations were applied to meet variance homoscedasticity when required. Besides, a principal component analysis (PCA) was conducted to reduce the dimensionality of soil nutrient data into two linear axes explaining the maximum amount of variance. Because our focal trees were spatially arranged within fragments, we evaluated if they could be considered as independent replicates since our dataset could have a spatial autocorrelation structure. Before modeling, we checked for spatial autocorrelation in soil respiration and enzymatic activity (glucosidase, phosphatase and chitinase). To do so, we used Moran I function (ape library, Gittleman & Kot, 1990). We found no spatial autocorrelations in our response variables, and hence, we considered our sampling points to be independent replicates. After our analysis, we performed a Moran test on the ANOVA model residuals and found no significant effects. In the case of SEM models, the standardized residuals were always <2.58, indicating non-significant discrepancy among variables (Grace, 2006). Thus, model residuals supported our initial assumption of sample independence. Soil functional, biotic and abiotic environmental variables were analyzed by two-way Analysis of Variance (ANOVA) considering the factors matrix influence (low, medium, high) and coverage (under canopy, open areas). Since coverage had a predominant effect on soil functioning, we assessed coverage and matrix effects separately. To do so, we firstly quantified the effects of tree cover for a given matrix level, and, secondly, the effects of matrix influence for a given coverage. For this purpose, we used one-way ANOVA. Tukey’s HSD were used as post hoc test (p < 0.05).

Structural equation models (SEMs) were used to test the direct and indirect influence of both biotic and abiotic factors on measured soil microbial functioning variables R_s and enzymatic activity (as latent variable). The latent variables are those not directly observed (i.e., enzymatic activity) but rather inferred from other observed variables, directly measured, denominated indicators (i.e.,  β-glucosidase, chitinase, and phosphatase). Our models considered a complete set of hypotheses based on literature, previous exploratory analyses and our own experience (Flores-Rentería et al., 2015; Flores-Rentería et al., 2016). In short, we hypothesized that matrix influence would promote nutrient input and tree growth (Flores-Rentería et al., 2015; Morán-López et al., 2016). Increased tree size would result in higher organic matter available in the soil (SOM), which in turn would modify abiotic conditions (i.e., increased moisture) (Abu-Hamdeh, 2001). These new abiotic conditions combined with increased SOM would result in higher respiration rates and overall enzymatic activity (Pérez-Izquierdo et al., 2017; Schlesinger & Andrews, 2000). Several models were run, and the best-fitted ones were finally selected according to the covariance proximity between observed and expected data (goodness-of-fit χ²). Standardized path coefficients were estimated by using the maximum likelihood algorithm (Shipley, 2002). The degree of fit between observed and expected covariance structures was assessed by root mean square error of approximation statistic (RMSEA) (Steiger, 1990). RMSEA values <0.08 indicate a good fit, between 0.08 and 0.10 provide a moderate fit, and >0.10 suggest a poor fit (Maccallum, Browne & Sugawara, 1996). Model fit to data was additionally evaluated by the goodness-of-fit index (GFI) and the Bentler and Bonett’s normed-fit index (NFI), both with values ranging between 0 and 1, and those >0.9 indicating an acceptable fit (Iriondo, Albert & Escudero, 2003). All statistical analyses were performed by using SPSS® and SPSS® AMOS 20.0 software’s (IBM Corporation Software Group, Somers, NY, USA).

Soil characteristics

Soil biotic and abiotic variables were strongly influenced by the canopy cover (Table 1). Under canopy, significantly higher values of soil moisture, SOM, and lower values of soil temperature and pH, were found compared to open areas (Table 1; Table S1). Matrix effects were related to soil temperature and moisture. Soil temperature tended to be lower in forest edges (intermediate influence), both under canopy and in open areas. In the case of moisture, matrix effects were mediated by coverage. Only in open areas it was significantly affected by the agricultural matrix, with the highest values recorded at forest edges (Table 1; Table S1). Organic matter and pH were not significantly affected by matrix influence. In areas with higher matrix influence, trees were higher (basal area, height and canopy projection), which resulted in an overall higher tree influence (Fig. S1). According to the PCA of soil variables, the PC1 axis reflected a gradient of nutrient availability related with the degree of forest fragmentation, with higher amounts of C, N, P, Ca²⁺, S, and Mg²⁺ in small fragments and at forest edges than in forest interiors (Fig. S2). Regarding biotic variables, tree cover affected microbial biomass but did not affect community composition (bacterial and fungal richness). In contrast, matrix influence did not modify microbial biomass but affected bacterial richness, with more OTUs at small fragments and forest edges in comparison with the forest interiors (Table 1, Table S1).

Matrix and tree cover influence in soil functioning

Tree cover significantly and positively affected all enzymatic activities, independently of the level of influence of the matrix (Table S1; Figs. 1A–1C). Maximum β-glucosidase (14.62 pmol mg⁻¹ min⁻¹), chitinase (2.87 pmol mg⁻¹ min⁻¹) and phosphatase (12.57 pmol mg⁻¹ min⁻¹) activities were obtained under canopy in small fragments (Figs. 1A–1C). Regarding matrix influence, its effects depended on the particular enzymatic activity. β-glucosidase activity was significantly lower at forest interiors, both under canopy and in open areas (Fig. 1A). Nonetheless, matrix effects on the rest of enzymatic activities were modulated by tree cover. While in open areas chitinase and phosphatase activities were unaffected by matrix influence, under canopy, these enzymatic activities were significantly higher in small fragments (Figs. 1B–1C). Fragmentation effects on soil respiration patterns were weaker. The agricultural matrix did not have a significant effect on R_s in any case, while tree cover did at forest edges, where R_s was significantly higher in open areas than under canopy (Fig. 1D).

The structural-equation model proposed for soil functioning (Fig. 2) provided a good fit as indicated by the non-significant f value (χ² = 29.7; p = 0.99) and by the goodness-of-fit indices (RMSEA < 0.001; NFI and GFI > 0.96). The indicators of the proposed latent variable of enzymatic activity (β-glucosidase, chitinase, and phosphatase) showed a good fit. All tested soil enzymatic activities showed high maximum likelihood estimates: phosphatase 0.93, β-glucosidase 0.89, and chitinase 0.61. Squared multiple relations also exhibited high amounts of explained variance of the different variables included in the model, especially for the soil functioning indicators: β-glucosidase (89.8%), soil enzymatic activity (82%), chitinase (77.1%), phosphatase (60.1%) and R_s (17.7%) (Table 2).

The proposed general fitted model showed that the agricultural matrix influence had a positive impact on PC1 axis (i.e., soil nutrients), tree influence, and soil temperature (Fig. 2). Tree influence significantly and directly affected many different soil properties such as soil temperature (negatively), SOM amount, nutrients, pH and microbial biomass (positively) (Fig. 2), which modified soil respiration and enzymatic activities.

According to standardized regression weights, microbial biomass was directly affected by SOM (positive), soil moisture (positive) and pH (negative) (Fig. 2). Soil enzymatic activity was mainly driven by microbial biomass, followed by the tree influence (positive) and soil pH (negative) (Fig. 2; Table 3 direct effects). Specific enzymatic activities were regulated by different pathways. β-glucosidase activity was modulated by matrix influence and soil respiration. Chitinase activity was positively related to soil moisture but negatively related to fungal richness. In the case of phosphatase, it was negatively directly affected by tree influence though total effects of this variable turned positive (T = 0.24, Table 2). Contrary to our expectations, neither SOM amount or nutrient input directly affected overall enzymatic activity (Fig. 2). However, in the case of SOM some indirect paths were observed, i.e., through its effects on microbial biomass and soil pH.

Soil respiration was positively and directly affected by the environmental factors, temperature and moisture, which were in turn indirectly influenced by trees through increased SOM input (Fig. 2, Table 2). As expected, soil temperature and moisture were negatively correlated. R_s was neither affected by microbial biomass, bacterial and fungal richness, SOM or tree influence (Fig. 2).