Study site

The study was carried out in the Navahondona forest (Cazorla mountain range, SE Spain). The area forms part of an abrupt calcareous mountain with a Mediterranean-type climate of wet and cold winters, and hot and dry summers. Mean annual rainfall is 1070 mm and mean annual temperature is 12 °C (based on climatic data from Nava de San Pedro Meteorological Station located within the Navahondona forest, 1290 m.a.s.l., 37º 52’N – 2º 53’W). In Cazorla mountain range, the altitudinal gradient for mean annual temperature (the adiabatic lapse rate) decreases linearly by 0.64 °C for every 100 m of altitude, and annual precipitation increases by ~84 mm for every 100 m of altitude (Linares & Tíscar, 2010). Consequently, altitude represents a suitable proxy of climatic variables in the area.

The main tree species in the Navahondona forest are Pinus nigra Arn. salzmannii (Dunal) Franco, Pinus pinaster Aiton., Quercus ilex L. and Quercus faginea Lamark. Respectively, these species represent 62%, 12%, 19% and 4% of current tree density in the forest. They are all native species to the area and regenerate naturally. The most abundant shrub species belong to the families Lamiaceae and Leguminosae (Tíscar, 2004). The Navahondona forest has been regularly harvested over the 20^th century, following silvicultural prescriptions based on the selection method (Tíscar et al., 2011). It is currently included in Cazorla, Segura and Las Villas Natural Park, but timber extraction is still allowed.

Data set

A total of 298 plots were randomly selected from the 2010 Navahondona Forest Inventory database (2010NFI hereafter) to compile the study dataset. The selected plots were located over two opposing hills (study area hereafter), covering an altitudinal gradient of 739 m (from 973 to 1712 m.a.s.l., median = 1338 m.a.s.l.) and 5569 ha of forest (Figure 1). The 2010NFI consisted of 15m radius circular plots systematically distributed on a square grid of 200 m throughout the study area. From the 2010NFI database, the altitude (m.a.s.l.), aspect (8 categories in intervals of 50 grades, ^g), total tree cover (the proportion of the forest floor covered by the vertical projection of the tree crowns, %) and stand basal area per tree species (m² ha^–1) were noted for each 15m radius plot. Additionally, the composition of the shrub layer and the cover of each species (percentage of area occupied by the vertical projection of the shrub foliage, %), and tree regeneration (number of individuals of each species with height < 1.30 m, seedlings hereafter), all measured in a 5m radius concentric subplot within the inventory plot, were obtained from the 2010NFI database for each selected plot.

Figure 1. Geographical location of the Navahondona forest in Cazorla mountain range (SE Spain) and situation of the samplings plots (dots) across the elevation and aspect gradients. Contour lines are separated by 100 m. Numbers indicate altitude above see level of the nearest contour line.

Data analysis

First, aspect was categorized into two groups: South (126^g – 325^g) and North facing slopes (326^g – 399^g and 0^g – 125^g), which respectively took the values zero and one in the statistical analysis. Then, I computed total shrub cover (%), shrub species richness (total number of species, S) and shrub diversity (Shannon diversity index, H’) for each plot. The Shannon diversity index was calculated as

where pi is the percentage cover of the i species divided by the total cover of the plot, and S is the number of species (Magurran, 1988).

A general linear model was used to analyse the effects of aspect, altitude and stand basal area on the percentage of shrub cover after the arcsine transformation of the original cover data. Likewise, the effects of the same three predictor variables: aspect, altitude and stand basal area on the dependent variables shrub diversity and shrub species richness were analysed by means of two Generalized Linear Models (GLM), using respectively a gaussian (link = identity) and a negative binomial (link = log) error distribution. Repeated estimates of diversity from the same area, as with inventory plots, are often normally distributed (Magurran, 1988) and, therefore, they can be handled in GLMs by specifying a Gaussian error distribution. In contrast, species richness represents count data, i.e. numbers of shrub species within each sampled plot, and this type of data usually have a Poisson distribution that can be modelled in GLMs by specifying a negative binomial error distribution (Crawley, 2013). Given that a strong correlation existed between tree cover and stand basal area (Spearman rho = 0.64, p < 0.0001), the effect of tree cover was not tested in order to avoid multicollinearity. Tree cover was dropped and stand basal area was included as an explicative variable in all the analysis, because basal area is measured and generally used to describe forest structure, while tree cover is typically estimated by sight in forest inventories. Stand basal area was included as a linear and second-order polynomial term to select the best transformation of this variable to account for nonlinearity (Ruiz-Benito et al., 2012).

I also fitted two separate GLMs with a log link function and a negative binomial error distribution to examine differences in mean Pinus nigra spp. salzmannii (Pinus nigra hereafter) and Quercus seedling density across topographic and stand structural gradients. Hereafter, the term Quercus will refer to the oak species: Quercus ilex and Quercus faginea grouped together. Regeneration of other tree species was negligible and, for this reason, not analysed. The predictor variables included in the full models were altitude, aspect, stand basal area, percentage of shrub cover and presence of a seed source (yes or no) of the target taxon: Pinus nigra or Quercus. In this respect, I considered that trees of both taxa with diameters at breast height (d.b.h.) ≥ 12.5 cm had attained sexual maturity. Therefore, the presence of a seed source was noted as positive when at least one tree of that size was enumerated in the inventory plot. The continuous variables (stand basal area and shrub cover) were included as a linear and second-order polynomial term to select the best transformation of these variables to account for nonlinearity (Ruiz-Benito et al., 2012). Due to the presence of overdispersion, a negative binomial distribution was preferred over a Poisson distribution in all the GLMs that implied a response variable in the form of count data (Crawley, 2013). Finally, I computed response curves of seedling density to altitude and stand basal area from the linear predictors of the best GLMs. In other to do this, I fixed the values of the other continuous variables at the observed mean and the most common value for categorical variables. Then, I back-transformed the expected values of seedling density from the resulting model between the minimum and maximum values of altitude or stand basal area to produce graphs with seedling density on the original scale.

In all the GLM analysis, the minimal GLM model was obtained with a step-wise procedure using the Akaike Information Criterion (stepAIC library MASS in R). A Canonical Correspondence Analysis (CCA) was conducted with function cca in library vegan of R (Oksanen et al., 2010) to relate the specific distribution of the 10 most abundant shrub species in the study area to gradients of altitude, aspect and stand basal area. Difference in mean Pinus nigra and Quercus seedling abundances, which could be interpreted as an indirect indication of successional trends, were examined using a non-parametric Mann-Whitney U test. Throughout the paper, values are means ± standard errors.

Shrub composition and diversity

2010NFI surveyors did not always identified shrubs at the species level and, frequently, they only noted the genus corresponding to each shrub present in the inventory plot. A total of 31 taxa (either identified at the genus or species level) were listed in the reports of the 298 sampled plots. The two taxa with the highest frequency in plots were Rosa sp. and Thymus sp.; and the two taxa with the highest percent cover in the plots where they were present were Rosmarinus officinalis and Echinospartum boissieri (Table 2). The first two CCA axes explained over 94% of the variation observed in percent cover of the 10 most abundant shrub species. The CCA (Figure 2) revealed that the presence of Lavandula latifolia is associated with North-facing slopes and that Echiospartum boissieri, Berberis vulgaris and Erinacea anthyllis are species from sites located at higher altitude, whereas Rosmarinus officinalis and Cistus sp. are mostly found on South-facing slopes and lower altitude. Most taxa thrived in plots with intermediate mean values of stand basal area, although Echinospartum boissieri was found in relatively open sites and Rubus ulmifolius in the densest plots (Table 2).

Figure 2. Canonical Correspondence Analysis (CCA) of the topographic and stand-structural variables characterizing the sites where the 10 most abundant shrub taxa in the Navahondona forest are found.

The three variables used to describe shrub communities attained low values in the Navahondona forest. Mean percent shrub cover was 6.88 ± 0.57%, mean shrub richness was 2.82 ± 0.10 species, and the mean value of the Shannon diversity index was 0.76 ± 0.03 nits per individual. As a consequence, shrub cover, richness and diversity were generally low across habitats (Table 1). Stand basal area was the only independent variable that significantly explained the observed variability in all three shrub responses, i.e., percent cover, richness and diversity (Table 3). The parameter estimates of stand basal area were negative for the three response variables and the introduction of a second-order term never improved the model fits (Table 3), indicating that shrub cover, richness and diversity might decrease linearly along a gradient of increasing tree density in the study area. Yet shrub cover was generally low and highly variable in open sites; specifically, percent shrub cover ranged from 0% to 54% (mean = 7.6 ± 1.33%) in plots with stand basal area ≤ 10 m² ha^–1 (n = 55). Conversely, mean percent shrub cover was as low as 2.04 ± 0.35% in plots with stand basal area ≥ 39 m² ha^–1 (n = 23).

Patterns of tree regeneration

Both Pinus nigra and Quercus regeneration were positively affected by the presence of seed-trees (pine and oak individuals with d.b.h. ≥ 12.5 cm), although the effect was only marginally significant for Pinus nigra recruitment. Quercus seelings were more abundant than Pinus nigra seedlings (5.23 ± 3.87 and 0.96 ± 1.7 individuals per plot respectively; Mann-Whitney U test, U = 14192, n = 298, p < 0.0001). The topographic variables altitude and aspect explained both Pinus nigra and Quercus recruitment patterns (Table 4). The density of Quercus seedlings decreased with increasing altitude, while Pinus nigra recruitment was favored by an increase in altitude (Figure 3a). In turn, the density of Quercus seedlings was lower on north-facing slopes than on south-facing slopes, but Pinus nigra seedlings benefited from the environmental conditions found on north-facing slopes, as denoted by the signs of the regression parameter estimates (Table 4).

Figure 3. Form of the relationship between predicted Pinus nigra recruitment (continuous line) and Quercus recruitment (dashed line) along gradients of altitude (a) and stand basal area (b) . Pinus nigra equations are: y = exp(-3.27+0.002x) and y = exp(-1.38+0.165x-0.003x2) for altitude and stand basal area, respectively. Quercus equations are: y = exp(5.56-0.003x) and y = exp(0.912+0.035x-0.0004x2) for altitude and stand basal area, respectively.

With respect to stand basal area, this variable affected both Pinus nigra and Quercus regeneration (Table 4). The value of Pinus nigra stand basal area increased with altitude (r = 0.32, n = 298, p<0.001), although individuals of this species were found throughout the study area (252 out of 298 sampled plots). On the contrary, the stand basal area of Quercus decreased with altitude (r = -0.24, n = 298, p<0.001), and trees of this taxon were found in only 166 sampled plots. Total stand basal area did not vary with altitude (r = 0.06, n = 298, p>0.10).

A nonlinear response of Pinus nigra recruitment was observed along the gradient of stand basal area and, similarly, second-order polynomial transformations of stand basal area and shrub cover were included in the best GLM for Quercus regeneration (Table 4). Response curves of Pinus nigra and Quercus recruitment along the tree density gradient (the range of the variable stand basal area was 0 – 71.85 m² ha^–1, mean ±SE = 22.24 ± 0.77 m² ha^–1) suggested that recruitment might be enhanced at intermediate values of stand basal area, whereas high values of this explicative variable exert a negative effect on seedling density. Specifically, seedling density was expected to be maximized at 27.50 m² ha^–1 of stand basal area for Pinus nigra, and at 43.75 m² ha^–1 for Quercus (Figure 3).

In contrast, the parameter of the second-order term was positive for the variable shrub cover (only included in the best GLM for Quercus regeneration (Table 4)). A response curve (not shown) indicated that this result implies a higher rate of survival at the top extreme of the shrub cover gradient.

Shrub composition and diversity

The relationship between shrubs and topographic variables has not been widely investigated in Spain. A recent study reported a regional negative effect of altitude on shrub cover in Pinus nigra forests from NE Spain (Coll et al., 2011), and, before this, Gracia et al. (2007) found that shrub richness and diversity tend to decrease with altitude in forests dominated by Pinus sylvestris and Pinus uncinata in the Pyrenees. To my knowledge, no other research has investigated this topic. In the present study, shrub cover, richness and diversity showed no significant variation with altitude in a Mediterranean forest located in southern Spain (Table 3). This disagreement with previous studies could be explained by the facts that (i) the area studied in the Navahondona forest was relatively small and did not include a latitudinal gradient, (ii) although the altitudinal gradient was strong (739 m), the most abundant shrub species were all Mediterranean and, consequently, shrub richness and diversity scores were not positively affected by the presence of a mixture of species with disparate biogeographical origins, (iii) and, compared with forests from northern Spain, shrub cover was much less in the study area (Table 1). Indeed, Pinus nigra forests from southern Spain attain the least shrub cover within the whole distribution area of this pine species (Tíscar & Linares, 2013). This circumstance could be related to forest structure according to results from the present study, but seed dispersal and herbivory could also be involved.

Although shrub recruitment might be nil in Mediterranean forests during very dry years (Mendoza et al., 2009b), it is unlikely that climatic conditions and, specifically, water availability have limited the expansion of shrub communities in the Navahondona forest for at least three reasons. First, results showed that altitude (a proxy of climate) does not influence shrub abundance, i.e. percent shrub cover does not increase as lower temperatures reduce evaporative demands as altitude increases. Second, shrub species most frequently found in Navahondona are Mediterranean (Table 2) and, consequently, well adapted to the prevalent climate in the area. Moreover, many of them are pioneer species resistant to arid conditions. Finally, climatic conditions should be suitable to maintain a larger biomass of shrubs, since a dense forest thrives in the region. Indeed, seed dispersal, herbivory and competition rather than rainfall are thought to control shrub establishment in a semi-arid area from SE Spain, which receives much less annual precipitation than the Navahondona forest (230 mm vs. 1070 mm (Pugnaire et al., 2006). Consequently, I will discuss how stand basal area (i.e. competition with canopy trees), seed dispersal, and herbivory might explain the characteristics of the shrub communities in the study area.

Results from this study were in agreement with previous studies that reported the existence of one-sided competition for water, nutrients and light between canopy trees and shrubs in Mediterranean mountain forests (Coll et al., 2011; Gómez-Aparicio et al., 2009). For instance, Figure 2 indicated the existence of competition for light since pioneer species, such as Echinospartum boissieri and Erinacea anthyllis, occupy sites with low tree density (or stand basal area), whereas shade-tolerant late successional taxa, such as Rosa sp. and Helleborus foetidus, are found at intermediate levels of stand basal area.

However, regardless of the species, shrubs were very infrequent in the densest forest stands, probably because of both deep shade conditions and canopy tree competition for soil nutrients and water (Gómez-Aparicio et al., 2009). High tree densities can additionally act as an obstacle to seed dispersal by wind and deter birds, limiting seed dispersal and impoverishing shrub communities inside dense forest stands (De la Montaña et al., 2006).

Poor seed dispersal, together with herbivory, could also explain scant shrub cover in sites with low values of stand basal area, where competition with canopy trees is reduced. The species listed in Table 2 have either fleshy- or dry-fruits. Seeds of the fleshy-fruited taxa are mainly dispersed by birds under plants of the same species (Herrera, 1984a), whereas dry-fruited shrubs disperse seeds within short distances (Herrera, 2000). As a result, the capacity of these shrubs to find and colonize new suitable habitats is rather limited. Indeed, new cohorts of pioneer shrubs establish in forest gaps located within the study area when seeds are directly sown onto the ground (Herrera, 2002); i.e., seed sowing at unoccupied sites resulted in recruitment, this being an indication of limited seed dispersal (Ehrlén & Eriksson, 2000). Herbivory should also be considered as a limiting factor for shrub expansion, since the Cazorla mountain range has typically been overgrazed by domestic and wild ungulates (García-González & Cuartas, 1989; Soriguer et al., 2003). In fact, although large herbivores usually prefer eating other shrubs to species belonging to the genuses Thymus, Lavandula or Rosmarinus (three of the most frequent shrubs in the study area), plants of these taxa are still consumed by the aforementioned animals (Fernández-Olalla et al.,2006).

The Navahondona forest along with the Cazorla mountain range harbour a very rich vascular flora (Gómez-Mercado, 2011). In this sense, the 2010NFI database listed no less than 31 taxa of shrubs with disparate ecological requirements, but shrub richness and diversity scores were generally low in the study area (Table 1). This contradictory result could be explained by the strong negative effect of stand basal area on shrub cover and by the association between diverse topography and niche separation of shrub species. Thus, changes in shrub composition along topographic gradients resulted in a remarkable beta-diversity, as shown by Figure 2, but alpha-diversity remained low at the plot level (Table 1).

Patterns of tree regeneration

Results showed that Pinus nigra regeneration is more abundant on north-facing slopes and at higher altitudes. Sites with these topographical characteristics typically represent the coolest and rainiest habitats within mountainous areas; additionally, the lower temperatures and direct radiation results in lower potential evapotranspiration (Linares & Tíscar, 2011). Previous studies have found that summer drought is the main cause impeding Pinus nigra regeneration at the local (Tíscar & Linares, 2011) and regional scale (Tíscar & Linares 2013). Therefore, Pinus nigra seedlings may establish more successfully on shady aspect and higher altitude, where levels of soil-water availability are expected to be higher.

In contrast to Pinus nigra, altitude exerted a remarkably high negative effect on Quercus regeneration (Figure 3). An immediate explanation for this result is that oak juveniles might be wounded by freezing temperatures and, particularly, late-winter frosts (Valladares et al., 2008). Temperatures typically descend with increased altitude (Linares & Tíscar, 2010). However, the observed patterns of Quercus regeneration could also be a reflection of current oak distribution in the study area, since Quercus populations rarely thrive above 1400 m.a.s.l. in the Navahondona forest, partly due to the objectives of forest management during the first half of the 20^th century (Tíscar, 2014). These findings and differences in seedling abundance (Table 1) suggest that currently dominant Pinus nigra forests could be replaced by more mixed pine-oak stands at lower altitude sites, whereas it would probably persist as the sole or dominant species at higher altitudes. A replacement of pine by oak would agree both with previous studies conducted in Mediterranean pine forests (Lookingbill & Zavala, 2000) and with the observed higher shade tolerance of Quercus compared with pine seedlings (Figure 3). However, successional change is likely to depend on the silvicultural treatments applied and on acorn dispersal. Thus, active management would be required to favor oaks over pines where pine is currently abundant, because Pinus nigra grows taller and lives longer than Quercus species (Ruiz de la Torre, 1979).

On the other hand, results from this study showed a positive relationship between Quercus regeneration and the presence of a seed source. This finding was interesting, because acorns can be dispersed over long distances by the European jay, Garrulus glandarius (Gómez, 2003). Unfortunately, there are no density estimates of jay populations for the study area, but this bird definitely inhabits the Navahondona forest (Martí & del Moral, 2004). Results consequently suggest that the role of jays in oak regeneration is somehow diminished in the study area compared to other sites from the same natural region (Gómez, 2003), since the establishment of oak seedlings is more likely immediately adjacent to conspecific adult trees in the region studied. Other animals, such as rodents, help oak regeneration by caching acorns in the soil, but dispersion is carried out to distances that are well within the radius of the study circular plots (Gómez et al., 2007). Therefore, their role as possible drivers of forest succession is probably limited in the Navahondona forest. Pinus nigra produces light, winged seeds that can be blown long distance in order to colonize new sites and, thus, the relationship between seedling density and the presence of a seed source within the sampled plot was only marginally significant (Table 4).

Recruitment of both Pinus nigra and Quercus was largely affected by tree density. The observed correlation between seedling density and stand basal area exhibited a hump-shaped pattern in both taxa (Figure 3b), indicating the existence of mechanisms associated with resource availability and stress conditions along the basal area gradient. Specifically, seedlings may benefit from the shady conditions beneath tree canopies, because shade lowers evaporative demands (Valladares et al., 2008) and, where measured, soil-moisture content has proved to be higher beneath pine trees than in open areas during the summer, the most critical season for seedling survival (Castro et al. 2004, 2005). However, the potential benefits of the pine overstory can be quickly reversed along a gradient of increasing stand density by the negative effect of competition for water (Gómez-Aparicio et al., 2009). As a result, recruitment of Pinus nigra and Quercus diminished at both ends of the tree density gradient (Figure 3b). The density of Quercus seedlings peaked at a higher value of stand basal area, probably, because this species exhibits higher shade tolerance than pines.

Statistical analysis did not detect a facilitative effect of shrubs on Pinus nigra regeneration (Table 4). Tíscar & Linares (2011) observed that Juniperus communis plants facilitate Pinus nigra recruitment in a location next to the study area. However, facilitative effects of shrubs on tree seedlings tend to be species-specific (Gómez-Aparicio et al., 2004), and the results reported here are congruent with the week nurse function observed for the most abundant shrub species in the study area (Siles et al., 2008). In contrast, shrub cover benefited Quercus regeneration, and this effect increased along the shrub cover gradient as indicated by a positive quadratic term in the polynomial transformation of the variable (Table 4) and a response curve (not shown). A plausible explanation for this finding is that ungulates avoid grazing patches of unpalatable shrubs, and the larger the patch the more noticeable the avoidance (Baraza et al., 2006). Otherwise, wild and domestic ungulates would typically feed on Quercus seedlings (García-González & Cuartas, 1989). In fact, shrubs with highest occupancy in the study area are unpalatable (Fernández-Olalla et al.,2006).

Implications for forest management

The most abundant shrubs identified by surveyors of the 2010NFI data were pioneer plants inhabiting degraded land. As a management option, these sites could be restored to late seral forests, unless these shrub species are considered part of the ecosystem and their presence worthy as contributors to important forest ecosystem processes. For instance, species of the Lamiaceae family are visited by rich assemblages of pollinating insects (Herrera, 1989) and, therefore, contribute remarkably to forest biodiversity. Similarly, fleshy-fruited shrubs are important for avian communities, since two-thirds of the passerine species present in Cazorla mountain range feed on shrub fleshy-fruits over the winter, with fruits of Rosa and Berberis eaten by birds of the genus Turdus in particular (Herrera 1984a, 1984b). Although the present study did not aim to analyse the demography of shrub species, results suggest the existence of limited recruitment that could be enhanced by a proper management of vegetation structure (Zamora et al., 2010). This recommendation of proactive management also applies to oak regeneration, if Navahondona is to become a more mixed forest in accordance with the currently established management plans for Cazorla, Segura and Las Villas Natural Park. Nevertheless, regeneration levels of Pinus nigra are even lower than those of Quercus species (177 ± 34 seedling ha^–1 as best (Table 1)). Such a figure of seedling density is considered insufficient to restock managed pinewoods (Serrada et al., 2008) and, consequently, might be indicating regeneration failure. However, no regeneration silvicultural treatments have been implemented in the forest during the last 15 years. It is not possible to isolate from the current database the relative contribution of individual factors that are known to influence Pinus nigra recruitment. However, results confirmed previous general knowledge about the effects of drought and stand basal area on seedling establishment (Tíscar & Linares, 2011). The correct management of stand structure in Pinus nigra stands appears to be more important than ever considering biodiversity implications in the context of changing climates.