Tree patch distribution drives spatial heterogeneity of soil traits in cork oak woodlands

Abstract

Key message

Spatial heterogeneity of soil resources is linked to tree patch distribution in Mediterranean cork oak woodlands. Tree patch distribution modifies soil traits by varying litterfall inputs due to different tree covers.

Context

The spatial heterogeneity of soil resources affects the structure and functioning of the different plant communities on ecosystems. In Mediterranean oak woodlands, the scattered trees play a key role in this spatial heterogeneity and might strongly influence ecosystem functioning and its productivity.

Aims

To assess the influence of the spatial pattern of trees and litter nutrient content on the spatial heterogeneity of soil properties.

Methods

We used a combination of geostatistical techniques and a linear mixed model to evaluate the spatial heterogeneity of soil and the seasonal and spatial variability of litter nutrient content, respectively.

Results

Soil parameters showed a high spatial heterogeneity. Tree cover was positively related with soil pH, and the organic matter, N, K, and Ca soil content. The return of nutrients to the soil via leaf fall had a marked seasonality and a high spatial variability, but this spatial variability had no effect on the spatial pattern of soil resources.

Conclusion

The spatial heterogeneity of soil in cork oak woodlands is mainly driven by tree patches distribution. The importance of the spatial heterogeneity of soil resources and the spatial pattern of trees on the functioning of the dehesa ecosystem makes it necessary to include them in plant nutrition studies and modeling approaches in these ecosystems.

1 Introduction

Cork oak (Quercus suber L.) is a sclerophyllous evergreen oak and one of the most important forest tree species of the Mediterranean basin. In southwestern Europe and northwestern Africa, most of oak woodlands have been traditionally managed for centuries as savannah-type ecosystems, called dehesas in Spain and montados in Portugal. Cork oak savannas support high levels of biodiversity and ecosystem services, and they represent a model of a sustainable ecosystem management with coexisting human activities and natural resource conservation (Bugalho et al. 2011). This led to its inclusion in the Annex I of the European Union Habitats Directive (92/43/CEE). The most valuable product of cork oak savannas is cork that is harvested at 9 to 14 years intervals, but other uses such as livestock raising, firewood, and agricultural crops coexist in these ecosystems (Olea and San Miguel-Ayanz 2006).

The spatial heterogeneity of soil resources has an impact not only on plant individual functioning (Quilchano et al. 2008) but also on the structure and functioning of plant communities and ecosystems (Ettema and Wardle 2002; Sardans and Peñuelas 2013). Oak savannas have been established through the thinning of the natural oak forest, and they are characterized by a sparse tree cover and a diversity of understory vegetation that ranges from shrubs to grasses (Olea and San Miguel-Ayanz 2006). Under this structure, scattered oaks have a positive effect on chemical and physical soil fertility by the accumulation of organic matter (Eviner and Chapin 2003; Gallardo 2003; Moreno and Obrador 2007; Moreno et al. 2007; Rhoades 1997) and improving soil moisture, water availability and soil enzyme activity (Gallardo et al. 2000; Garcia et al. 2002). This spatial pattern allows the coexistence of plant species with different competition abilities (Ludwig et al. 2004; Moreno et al. 2005, 2007) and ultimately improves the production, quality, and diversity of herbaceous species (Marañón 1986). However, most of the studies concerning the impact of tree spatial distribution on soil heterogeneity have been carried out in holm oak savannas with a tree density of 10–50 trees per hectare, which are mainly used for livestock raising and crop production. Less attention has been paid to less disturbed and denser systems (over 100 trees per hectare) where soil features are expected to be more homogenized than in sparser woodlands (Lister et al. 2000; Quilchano et al. 2008). Increasing our knowledge about the spatial heterogeneity of soil resources in these systems is, therefore, of great interest for the selection of new areas with the potentiality of being intercropped. Furthermore, most of the studies about the spatial variability of soil resources in these ecosystems have been focused in the uppermost soil layer. Surprisingly, there is a lack of information about this variability at deep layers, despite the fact that oaks growing in dehesas show a deep root system with a low root density in the upper 10 cm of the soil (Moreno et al. 2005). In addition, there is also limited information about the spatial variability of litter nutrient content and its relationship with the spatial distribution of soil properties. Soil resources may determine the mineral content of leaves and, therefore, the characteristics of leaf fall, triggering a positive feedback between trees and soil (Aponte et al. 2011). Thus, increasing the knowledge of the spatial pattern of plant resources and processes affecting plant-soil interactions may contribute to a better understanding of the ongoing processes and subsequently have the potential to improve both management and conservation attempts.

In order to better understand tree-soil relations in Mediterranean managed cork oak woodlands, we aimed to (i) evaluate the spatial heterogeneity of soil chemical and physical properties at two soil layers (0–20 and 20–100 cm), (ii) to assess the effect of the spatial pattern of trees on these soil properties, (iii) to examine the seasonal and spatial variation of nutrient concentration in leaf fall, and (iv) to evaluate the effect of the spatial variability of litter nutrient content on soil properties. We hypothesized that (i) the spatial pattern of soil properties is mainly driven by tree cover in denser cork oak woodlands (with lack of cropping and/or livestock raising) and (ii) the tree cover influences soil chemical characteristics by the litter inputs.

2 Materials and methods

2.1 Site description

The study site is located in a Mediterranean cork oak forest close to the Doñana National Park (Huelva, SW Spain). We laid out a plot of 1.9 ha in a flat area at an altitude of 100 m a.s.l (37° 19′ 13″ N, 6° 25′ 37″ W) and with a tree density of 99.6 trees ha⁻¹. Cork oak is the dominant tree species, with a density of 83.5 trees ha⁻¹ and a basal area of 7.2 m² ha⁻¹, while holm oak (Quercus ilex subsp. ballota [Desf.] Samp) appears as secondary tree species (16.1 trees ha⁻¹ and 0.9 m² ha⁻¹) (Fig. 1). The understory community is mainly composed of Pistacia lentiscus, Daphne gnidium, Chamaerops humilis, Halimium halimifolium, and several Cistus spp. The main use of this forest is cork extraction, with no agricultural exploitation or livestock grazing. The understory management in the area consists in periodical and selective shrub clearing for fire prevention. The climate is typical Mediterranean with a mean annual precipitation of 579 mm, a mean annual temperature of 18.9 °C, a summer drought period of 5 months, and high interannual variability in temperature and precipitation (Pilas “Medina Garvey” weather station, at 11,300 m away from the study plot, data for the period 1966–2005). During the period studied (March 2004–February 2006), annual precipitation was 475 mm (153) from March 2004–February 2005 and 314 mm (67.8) from March 2005–February 2006 (spring precipitation in parenthesis). Soil was characterized by opening two soil profiles (online resource 1). The soil has a complex profile, with a sandy-loam to loamy-sand upper layer (top 25–40 cm) over an argillic horizon (with loam-clay to clay texture) that has stagnic properties and presence of irregularly distributed nodules of limestone. The soil is classified as Haplic Regosol over a Stagnic Regosol (IUSS 2007).

Fig. 1

Map of the study plot showing cork oak sampled trees (black dots), stripped cork oaks (open dots), non-stripped cork oaks (dots with cross), and holm oaks (framed dots) positions in the plot. Circle size is proportional to the DBH of the trees. Numbers represent the DBH (cm) for the sampled cork oaks. Soil sampling points are indicated with crosses

2.2 Field measurements and sampling

The positions of trees were established in a local coordinate system using a topographical total station. For each tree, tree diameter at breast height (DBH) and four crown radii, in the direction of each cardinal point, were measured. The range of DBH was 15.3 to 54.1 cm, the mean DBH was 28.3 ± 0.8 cm, and the mean crown diameter was 7.8 ± 0.2 m.

Litterfall samples were collected at 12 selected cork oak trees. For the selection, all trees of the plot were grouped in sets of three trees with similar diameters. Four of these groups, from those representing the central circumference class, were randomly selected to choose the 12 sampling trees. The DBH of sampling trees was 22 to 35 cm, and the mean DBH was 27.7 cm. Four 0.16-m² circular traps were placed at each cardinal point of every sampling tree at a distance corresponding to three quarters of the crown radius measured from the stem. Samples were collected monthly for 2 years in order to account for the temporal variability in nutrient return by litterfall before testing its relationship with the spatial distribution of soil parameters.

Soil samples were collected in the interception points of a 20 × 20-m regular grid established within the limits of the plot (Fig. 1). The sampling protocol resulted in 38 points that were sampled at the end of the litterfall sampling period (March 2006). Soil samples were collected with a drill every 20-cm depth until 100-cm depth. At each point, two subsamples were obtained: one superficial sample (0–20 cm) and one composed sample (20–100 cm).

2.3 Laboratory analysis

Litterfall samples were oven dried at 65 °C for 48 h, and cork oak leaves were separated from the rest of the plant material and weighed (±0.01 g). Leaves of each tree were ground with a micro mill to pass a 250-μm mesh. N was measured with an elemental analyzer (Thermo Finnigan 1112 Series EA, Milan, Italy). P, K, Ca, and Mg were measured with an ICP analyzer (ICP-OES, Jobin Yvon Ultima 2, Tokyo, Japan) after wet digestion with HCl (Temminghoff and Houba 2004).

Soil samples were air-dried (<30 °C) and sieved to pass a 2-mm mesh. For soil nutrient analysis, the Kjeldahl method was used for N determination, UV-visible spectrophotometry with Mehlich 3 extraction was used for P determination, and exchangeable K, Ca, Mg, and Na were determined by atomic absorption spectrophotometry with ammonium acetate extraction. In addition, soil characteristics were assessed according to recommended analytical procedures (IUSS 2007; Perkins et al. 2013). The parameters analyzed were soil particle size distribution by the pipette method, soil pH by a pH-meter (soil/water 1:5), soil organic matter (SOM) by the loss on ignition method, and cation exchange capacity (CEC) by the ammonium acetate method.

2.4 Data analysis

2.4.1 Spatial heterogeneity of soil properties

The spatial dependence of soil parameters was evaluated using geostatistical analyses. An empirical semivariogram was calculated for each parameter considering a lag size of 20 m (equivalent to the sampling grid side) and seven lags. The semivariogram was set to a spherical type function and then range, partial sill (C), and nugget (Co) were estimated. The spatial dependence of soil parameters was assessed by calculating the portion of the total variance (C + Co) explained by the spatial autocorrelation (i.e., C/C + Co).

Soil parameters were then interpolated to the whole plot by using ordinary kriging, and considering eight neighbors and only one sector. Prediction was isotropic and predominant directions were not observed. All geostatistical analyses were performed with ArcGIS ver. 9.2 (ESRI).

2.4.2 Relationship between tree cover and soil parameters

The effect of the presence of trees on soil properties is roughly proportional to their crown area projected to the soil surface (Zinke 1962). Hence, we calculated the total crown projected area (in m²) within a 10-m buffer at each soil sampling point in order to assess the effect of trees on soil parameters. Tree crown projected area was estimated for all trees as an ellipse calculated from the four measured crown radii (Sect. 2.2). The software ArcGIS ver. 9.2 (ESRI) was employed for this process. A buffer of 10-m radius was selected to avoid buffer overlapping, because the distance between adjacent sampling points is 20 m. Furthermore, the plot density is 99.6 trees ha⁻¹, so the mean distance between adjacent trees in the plot is around 10 m. The relationship between tree cover and soil properties was explored by regression analyses in which the total crown projected area was considered as explanatory variable. Normality of data was checked through the Kolmogorov-Smirnov test, and data was log-transformed when necessary. Statistical analyses were performed with SAS/ETS (ver. 9.2).

2.4.3 Model of nutrient concentration in leaf fall

We determined the nutrient content of leaf fall in 12 trees for each month over 2 years, so our data structure resulted in repeated measurements on each tree, and within-tree observations were autocorrelated. Hence, the hypothesis of the independence of observations that leads to a general linear model structure was not logical. Thus, we used a linear mixed model in which tree was considered as a random effect to assess the variability among trees in leaf fall nutrient content and to account for the correlation of within-tree observations. The seasonality of nutrient concentration in leaf fall was described as the month effect, the between-year variation as the year effect, and the seasonality between years as the interaction between month and year. Nutrient concentrations were transformed using the arcsin-square-root transformation. We also analyzed the seasonality in nutrient return among trees by adding the random interactions of tree × month and tree × year. Based on these considerations, the nutrient return to soil by leaf fall was evaluated using a linear mixed model with the following initial structure:

$$ {y}_{ijk}=\mu +{t}_i+{\alpha}_j+{\beta}_k+{\left(\alpha \beta \right)}_{jk}+{\left(t\alpha \right)}_{ij}+{\left(t\beta \right)}_{ik}+{e}_{ijk} $$

(1)

where y _ijk is the $\arcsin \sqrt{nutrient\left(\%\right)/100}$ transformation of nutrient concentration (N,P,K…in %) in fallen leaves from tree i of month k and year j, μ is the general mean, t _i is the tree random effect with i = 1,2,…,12, with the initial hypothesis t _i ~N(0, σ _t ²), α _j is the year fixed effect with j = 1,2 (year 1: March 2004 to February 2005, year 2: March 2005 to February 2006), β _k is the month fixed effect with k = 1,2,…,12, (αβ) _jk is the fixed interaction between the year and month, tα _ij is the random interaction between tree and year with tα _ij ~(0, σ _ij ²), tβ _ik is the random interaction between tree and month with tβ _ik ~(0, σ _ik ²), and e _ijk is the random error term under the hypothesis that e _ijk ~N(0,R).

We tested different variance-covariance structures of within-tree observations: autoregressive order 1, autoregressive heterogeneous, Toeplitz up to six bands, homogeneous Toeplitz up to three bands, unstructured up to six bands, and compound symmetry heterogeneous (see Littell et al. [2006] for a detailed description). Some of these structures are heterogeneous (i.e., variance among months is not homogeneous), because we could expect that months with higher values of the dependent variable have higher variability. Variance components for each structure were estimated by restricted maximum likelihood (REML) and the value of the statistic −2 × log likelihood (−2LL) was calculated. The Akaike information criterion (AIC, Akaike 1974), in which a lower value indicates a better model, was used for model selection. Models including and excluding the tree random effect and random interactions were tested by calculation of AIC. If the tree random effect was significant, we assessed the presence of spatial autocorrelation by considering an isotropic power covariance form. Finally, fixed effects were estimated by using generalized least squares (GLS), and the significance of each effect was determined with an F test. Only significant effects (α = 0.05) were retained in the model. The Scheffé test was used to compare significant effects.

The effect of nutrient inputs by litter on soil properties was analyzed by introducing soil variables as linear covariates in the models. For assignment of a particular value of the soil variable to each oak tree, we considered the area of influence of each individual to correspond to the Voronoi polygon in which the tree was included. An F test was used to assess the significance of soil covariates. A regression analysis was performed to explore the relationships between the significant soil covariates and the litter nutrient content. All statistical analyses were performed with SAS/ETS (ver. 9.2).

3 Results

3.1 Spatial heterogeneity of soil properties

The preliminary data exploration indicated large data variability for most of the soil parameters (Table 1). The coefficient of variation (CV) ranged from 0.10 to 0.92, with the highest CV values for the superficial CEC and exchangeable Ca of both soil layers, and lowest values for the pH and the percentage of sand of both soil layers (Table 1). Generally, the CV of each soil parameter was similar between soil layers or higher for the deep layer.

Table 1 Soil characteristics of the study plot for the upper 20 cm (n = 38) and for the layer 20–100 cm (n = 38)

The empirical semivariograms were successfully fitted to a spherical model for most of the parameters, showing spatial dependence within the plot (Table 2). Only the silt percentage of the shallow layer and the exchangeable K of both layers showed a random distribution within the plot (i.e., nugget model with C / (C + Co) = 0). The spatial dependence ranged from 0.10 to 0.88 (Table 2) with the pH, SOM, and exchangeable Ca of the superficial layer and the sand percentage, CEC and exchangeable Ca of the deep layer showing higher spatial dependence (C / (C + Co) > 0.6). Spatial ranges also highly varied among soil parameters (Table 2), but analogously to the CV, the spatial range of the parameters was similar between soil layers or slightly higher for the deep layer. Figure 2 shows the SOM, pH, N, and exchangeable Ca of the two soil layers as predicted by kriging. The predicted maps for the rest of soil parameters are shown in the online resource 2.

Table 2 Parameters of the spherical semivariogram models of soil properties

Fig. 2

Maps of organic matter content (%), pH, nitrogen content (g kg⁻¹), and exchangeable calcium (cmol⁺ kg⁻¹) predicted by kriging in the study plot. Panels in the left and right columns show prediction for the upper (0–20 cm) and lower (20–100 cm) soil layers, respectively. Polygons depict the projected crown area of all trees in the study plot

3.2 Relationship between tree cover and soil properties and nutrients

Total crown projected area was positively related with the SOM content (p = 0.009, R ² = 0.188, r = 0.433), pH (p = 0.009, R ² = 0.188, r = 0.434) and N content (p = 0.010, R ² = 0.189, r = 0.435) of the top 20 cm of soil. The exchangeable Ca of the top 20 cm of soil was also related with the crown projected area but at p < 0.10 (p = 0.067, R ² = 0.095, r = 0.308). In addition, crown projected area was also positively related with the N content and the exchangeable K of the 20–100-cm soil layer (p = 0.014, R ² = 0.165, r = 0.406; and p = 0.036, R ² = 0.123, r = 0.351, respectively).

3.3 Nutrient variability in leaf fall

The mean leaf fall (±standard deviation) for both years was 233.02 ± 158.57 g m⁻² year⁻¹, with greater leaf fall in year 1 than in year 2 (274.39 ± 166.56 g m⁻² and 191.65 ± 150.58 g m⁻², respectively). Leaf fall was greatest between March and May with a peak in April (Fig. 3, bottom right). The mean values (±standard deviation) of nutrient concentration in leaf fall were 8.69 ± 2.13 mg g⁻¹ for N, 0.32 ± 0.11 mg g⁻¹ for P, 4.79 ± 1.86 mg g⁻¹ for K, 8.10 ± 1.85 mg g⁻¹ for Ca, and 1.09 ± 0.27 mg g⁻¹ for Mg.

Fig. 3

Mean (±standard error) nutrient concentrations in leaf fall (mg g⁻¹) and leaf fall collected (g m⁻²) in each month of year 1 and year 2 (N = 12 trees)

Our initial analysis indicated that non-homogeneous variance structures performed better than homogeneous structures and that the introduction of the tree random effect significantly reduced the AIC and the −2LL (online resource 3). The percentage of variability accounted for by the tree random effect varied among nutrients and among months and was 3–35 % for N, 8–46 % for P, 31–60 % for K, 11–35 % for Ca, and 51–81 % for Mg. The interaction of tree × month was only significant for P, although the percent variance accounted for by this interaction was small (1–5 %). By contrast, the percent variance accounted for by the interaction of tree × year was higher for nutrients in which the interaction was significant: K (19–37 %), Ca (20–63 %), and Mg (12–19 %). In addition, the spatial covariance structure was significant for P, K, and Mg (online resource 3).

The interaction of year × month was highly significant (p < 0.0001, Table 3) for all nutrients. Figure 3 shows the nutrient concentrations in leaf fall at the year × month level. N and P have a similar seasonality, with lower concentrations during the peak of maximum litterfall (April) and increasing concentrations from the end of the summer until the next peak of litterfall. K also showed lower concentrations in April, but unlike N and P, K showed higher concentration during the summer. On the other hand, Ca showed the opposite trend with higher values during the peak of litterfall, whereas Mg did not show a marked seasonality like the other nutrients.

Table 3 Significance of the fixed effects for each selected model of nutrient concentration in leaf fall

In spite of the seasonality of nutrients in leaf fall that is rather similar between years, there were some significant differences for some nutrients (Fig. 3). N was higher in year 1 from June to August and from November to February. P also varied between years with generally higher concentrations in year 2, especially after the summer, and Ca also showed higher concentrations in year 2 from June to November.

3.4 Effect of litter nutrient content on soil properties

Soil covariates were only significant in the case of Mg litter content. The significant covariates were exchangeable Ca of the top 20 cm of soil (Ca_20) (p = 0.0067), exchangeable Ca of the layer 20–100 cm of soil (Ca_100) (p = 0.0003), and pH of the layer 20–100 cm of soil (pH_100) (p = 0.0353). The introduction of the covariate Ca_20 in the Mg model reduced the variance component at the tree level from 1545.16 to 361.99, so it accounts for 76 % of the variance at this level and for 39–61 % of the total variance, depending of the month. The covariates Ca_100 and pH_100 accounted for 54 and 68 % of the variance at the tree level, respectively. The estimated coefficient of the significant covariates and the regression coefficients were negative. In other words, lower Mg litter content was associated with higher soil exchangeable Ca and pH (Fig. 4), in part because soil Ca content and pH are correlated.

Fig. 4

Relationship of leaf fall magnesium concentration with soil exchangeable calcium at a depth of 0–20 cm (A) and with soil pH at a depth of 20–100 cm. The correlation coefficients (r), the p value, the coefficient of determination (R ²), and the equation for the regression lines are indicated

4 Discussion

4.1 Spatial heterogeneity of soil properties and its relation with tree cover

Trees are the main drivers of the spatial heterogeneity of soil in Mediterranean oak savannas, where soil patches strongly influence the spatial patterns of understory vegetation, community structure, and ecosystem function (Gallardo 2003; Moreno and Obrador 2007; Moreno et al. 2007). Our results confirmed the high spatial heterogeneity of soil properties in seminatural cork oak woodlands without livestock raising and/or crop production. We found that the spatial range (i.e., distance at which samples remain spatially dependent) was similar for the SOM, N content, and exchangeable Ca. The accumulation of SOM through litterfall is related with higher N mineralization (Gallardo 2003). In addition, the cork oak leaf fall is especially rich in Ca (Andivia et al. 2010); thus, the distribution of SOM, N, and Ca in soils is expected to be related. The spatial range of these parameters (85–120 m) roughly matched with the distance between the northern (higher tree density) and the southern (lower tree density) part of the study plot (Figs. 1 and 2), indicating that the spatial pattern of these parameters is related with the spatial pattern of trees. Indeed, we found a positive relationship between the tree cover and the SOM, N content, and exchangeable Ca in the soil that corroborates our initial hypothesis. We also found a positive relationship between the tree cover and the exchangeable K in the soil that could be related with the throughfall and K leaching from tree canopy (Gallardo 2003). However, we did not find a spatial dependence of this nutrient within our plot probably because the spatial autocorrelation of this element occurs at a lower spatial scale than the minimum scale sampled (20 m). Superficial pH and CEC showed similar spatial ranges but lower than SOM or exchangeable Ca which are expected to affect the spatial distribution of these parameters. Nevertheless, these parameters are highly dependent upon soil texture (Binkley and Fisher 2013). In fact, the spatial autocorrelation of superficial pH and CEC matched the spatial autocorrelation of superficial clay and sand content. Thus, at the study scale, the spatial pattern of pH and CEC seems to be more related to soil characteristics than to the spatial pattern of trees. However, the positive correlation between tree cover and pH might indicate that tree patch distribution has a greater effect on soil pH at lower scale. P and exchangeable Mg also showed similar range than CEC, pH, clay, and sand content. This confirms that P availability in soils depends not only on the balance between mineralization and uptake but also on the interaction with soil minerals, and probably on the distribution and density of roots and mycorrhizas (Gallardo 2003; Sardans and Peñuelas 2013). The likely effect of tree patches on soil P distribution might be hampered by the irregular distribution of shrubs in the plot and its effect on creating microenvironments (Rolo et al. 2012) and/or by the higher resorption efficiency of cork oak trees that leads to low P concentration in leaf fall (Andivia et al. 2010; Caritat et al. 1996; Robert et al. 1996). We did not find any relationship between tree cover and soil texture. The main effect of trees on soil structure is due to root disaggregation of rock fragments; thus, differences in belowground traits such as root morphology, C inputs, or associated soil biota are more likely to be related with the heterogeneity of soil texture than the mere presence of trees (Eviner and Chapin 2003).

Our results show the effect of tree cover on soil heterogeneity in denser cork oak woodlands. The spatial heterogeneity of soil resources is a key factor in determining the community structures of grasses in these ecosystems (Gea et al. 2009; Moreno et al. 2007). Modeling approaches considering the effect of tree patches on the spatial heterogeneity of soil and tree-crop competitive interactions might be of great interest to establish management guidelines for intercrops sustainability and profitability in these ecosystems.

4.2 Nutrient variability in leaf fall

The highly significant effects of month and the interaction of year × month (Table 3) indicate that the concentration of nutrients in leaf fall varies seasonally and is related to tree phenology (Andivia et al. 2010; Caritat et al. 1996; Robert et al. 1996). In general, Mediterranean plants reduce nutrient losses through litterfall by the retranslocation of nutrients from senescent leaves to new leaves before leaf shedding, with the exception of less mobile nutrients such as Ca (Aerts 1995; Escudero and Del Arco 1987; Oliveira et al. 1996). However, we found differences in the seasonality of leaf fall nutrient concentration between years, probably induced by the lower precipitations during the second year and the consequent higher leaf retention, as is corroborated by the smaller leaf fall peak in year 2 (Fig. 3). As a result, there was an accumulation of less mobile nutrients (e.g., Ca) and a reduction in the content of growth-limiting nutrients, such as N, in leaf fall during the second year, although this trend was not observed for P especially after the summer. The high seasonality of nutrient content in litterfall and its inter-annual variability highlight the need of considering this variability to more accurately estimate the annual influx of nutrients into the soil.

As a consequence of reducing nutrient losses through litterfall, litter is more recalcitrant, and C/nutrient ratios are higher leading to lower decomposition rates (Incerti et al. 2011; Sardans and Peñuelas 2013). This strategy enhances the control of nutrient stocks and the capacity of the tree-soil system to retain nutrients in the upper part of the forest floor that may become readily available to the trees under favorable soil temperature and moisture (Oliveira et al. 1996; Sardans and Peñuelas 2013). Nonetheless, the highest nitrification rates in soils of dehesas match the highest litterfall inputs (spring-summer), indicating its dependence on C availability (Gallardo et al. 2000) and therefore on the spatial and temporal distribution of SOM. In addition, seasonal changes on soil moisture have been identified as the main factor in controlling soil chemistry and activity in Mediterranean ecosystems (Sardans and Peñuelas 2013). Hence, cork oak nutrient uptake, nutrition, and ultimately, nutrient return to the soil through litterfall are strongly influenced by seasonality, and that temporal variation might represent a more determinant factor than its spatial variability (Orgeas et al. 2002). However, a large part of the variability in nutrient concentration of leaf fall has been explained at the tree level, and the spatial covariance was significant for P, K, and Mg.

4.3 Effect of litter nutrient content on soil properties

The positive relationship between tree cover and soil nutrient content together with the spatial autocorrelation of litter nutrient content might lead us to expect higher nutrient content in the soil close to those trees with higher litter nutrient content. Nonetheless, we did not find any significant effect of the spatial variability of litter nutrient content on soil properties. This might indicate that the variability in nutrient return by litterfall among trees is not enough to induce significant differences in soil properties beneath these trees and that the mere presence of trees induces the observed soil heterogeneity. However, we found that leaf fall Mg decreased with increasing soil Ca and soil pH (Fig. 4). Two hypotheses can be suggested. It is known that Mg²⁺ uptake is strongly depressed by high levels of soil Ca²⁺ due to competition for binding sites at the soil micelle surface and the rather low binding strength of Mg²⁺ (Marschner 1995). On the other hand, Mediterranean oaks have been proved to modify soil conditions leading to a positive tree-soil feedback by increasing the fitness of individuals of the same species (Aponte et al. 2011). Q. suber is a calcifugous tree that prefers acidic soils (Montero and Cañellas 1998); therefore, the lower leaf fall Mg content might be a strategy to cope with high levels of soil pH through decreasing the inputs of mobile cations to the soil. However, the low Mg concentration reported in the fresh leaves of these trees (Andivia et al. 2010) in comparison with other studies (Oliveira et al. 1996; Orgeas et al. 2002; Passarinho et al. 2006) suggests that the hypothesis of the Ca-Mg antagonistic effect is more realistic.

This antagonist effect and its negative consequences on the nutritional status of cork oak trees point out the importance of a thorough assessment of soil conditions to assist in the detection of potential forest species best adapted to site-specific conditions. In this direction, the establishment of new cork oak plantations in areas with high soil Ca content and pH might limit tree growth and fitness, and therefore the production of high-quality cork.

5 Conclusions

Soil parameters showed a high spatial heterogeneity in managed cork oak woodlands. Trees were the main source of spatial variability for those parameters related with the accumulation of organic matter through litterfall (i.e., SOM, N, and Ca content). On the other hand, the spatial pattern of pH and CEC seems to be more related to other soil properties, like soil texture, than to the spatial pattern of trees.

The variability of the nutrient return to the soil via leaf fall showed a high seasonal and spatial variability that should be taken into account when the annual influx of nutrients into the soil is calculated. However, the spatial variability of litter nutrient content did not induce differences in soil properties. Our results suggest that the spatial heterogeneity of soil resources and the spatial pattern of trees are key factors to be considered in plant nutrition studies and modeling approaches in these ecosystems.