| Citation: |
Álvaro Rubio-Cuadrado, Iciar Alberdi, Isabel Cañellas, Fernando Montes, Jesús Rodríguez-Calcerrada, Rosana López, Guillermo G. Gordaliza, María Valbuena-Carabaña, Nikos Nanos, Ramón Perea, Luis Gil. Twenty years of population dynamics in European beech-oak forest at their rear range margin anticipate changes in its structure and composition[J]. Forest Ecosystems, 2024, 11(1): 100197. DOI: 10.1016/j.fecs.2024.100197
|
There is an increasing interest in restoring degraded forests, which occupy half of the forest areas. Among the forms of restoration, passive restoration, which involves the elimination of degrading factors and the free evolution of natural dynamics by applying minimal or no management, is gaining attention. Natural dynamics is difficult to predict due to the influence of multiple interacting factors such as climatic and edaphic conditions, composition and abundance of species, and the successional character of these species. Here, we study the natural dynamics of a mixed forest located in central Spain, which maintained an open forest structure, due to intensive use, until grazing and cutting were banned in the 1960s. The most frequent woody species in this forest are Fagus sylvatica, Quercus petraea, Quercus pyrenaica, Ilex aquifolium, Sorbus aucuparia, Sorbus aria and Prunus avium, with contrasting shade and drought tolerance. These species are common in temperate European deciduous forest and are found here near their southern distribution limit, except for Q. pyrenaica. In order to analyze forest dynamics and composition, three inventories were carried out in 1994, 2005 and 2015. Our results show that, despite the Mediterranean influence, the natural dynamics of this forest has been mainly determined by different levels of shade tolerance. After the abandonment of grazing and cutting, Q. pyrenaica expanded rapidly due to its lower shade tolerance, whereas after canopy closure and forest densification, shade-tolerant species gained ground, particularly F. sylvatica, despite its lower drought and late-frost tolerance. If the current dynamics continue, F. sylvatica will overtake the rest of the species, which will be relegated to sites with shallow soils and steep slopes. Simultaneously, all the multi-centennial beech trees, which are undergoing a rapid mortality and decline process, will disappear.
Forests are estimated to cover 31% of the global land area, that is, 4.06 billion ha (FAO, 2020), of which nearly 2 billion ha are degraded (Minnemayer et al., 2011), needing major restoration efforts by governments and other institutions. Among the different types of restoration, passive (spontaneous) approaches, in which the only intervention is removing degrading activities to allow natural succession, are becoming increasingly common (Atkinson and Bonser, 2020). As long as there are nearby seed sources, the result of cessation of these activities generally leads to an increase in forest density until the canopy closes (Marchetti et al., 2018; Martín-Forés et al., 2020). The process of densification depends on many factors, such as the standing species distribution, their seed production, their shade tolerance, topography, climate, soil, etc. The complex interaction of these factors makes the evolution and dynamics of the forest difficult to predict.
One of the most important factors determining forest dynamics is shade tolerance. Colonizing shade-intolerant species will benefit in treeless or low-density areas. On the contrary, shade-tolerant species will benefit in areas with high competition for light and gap dynamics (McCarthy, 2001). In mixed forests, where species with different shade tolerances coexist, the dynamics in the absence of disturbance usually lead to the displacement of species with lower tolerance and the formation of nearly or completely monospecific forests, as commonly occurs in mixed beech forests (Petritan et al., 2017; Knapp and Pallardy, 2018; Rubio-Cuadrado et al., 2018). Another factor related to shade tolerance is the competitive ability of the species, i.e. the ability to maintain higher relative growth rates than other species under conditions of high competition. In this sense, beech also tends to have a higher competitive ability than the coexisting species (Hein and Dhôte, 2006; Manso et al., 2015; Rubio-Cuadrado et al., 2020a). The positive correlation between growth and vigor (Camarero et al., 2015) results in loss of vigor of species with less competitive capacity and, eventually their displacement by more competitive species.
An additional element that conditions the dynamics of densification is the adaptation of species to an ongoing rapid climate change. In southern Europe rising temperatures and increased drought, will negatively affect temperate species located at the southern limit (rear edge) of their distribution range (Hampe and Petit, 2005), favoring Mediterranean species better adapted to drier and warmer conditions (Rubio-Cuadrado et al., 2018; Perea et al., 2020). Other phenomena, like higher frequency of late frosts, can also alter the forest dynamics by shifting relative competitive abilities amongst the tree species (Sangüesa-Barreda et al., 2021; Rubio-Cuadrado et al., 2021a).
Another important factor to take into account is the regeneration capacity of the different species present prior to the cessation of the degrading activity. In this sense, the regeneration problems that some Quercus species tend to have are well known, whether due to lower seed production, high seed predation and herbivory, or competition with other plants in the understory (Shaw, 1968; Mellanby, 1968; Lorimer et al., 1994; Perea et al., 2020). This may result in these species not expanding or densifying where they are present (Marchetti et al., 2018). Soil characteristics (depth, nutrients, water retention capacity, etc.) also influence the dynamics of the forest, favoring the establishment and expansion of some species or others according to their requirements (Shibata and Nakashizuka, 1995; Masaki et al., 2005; Diaci et al., 2020). In degraded forests, loss of the fertile topsoil in eroded areas can hamper recruitment of species established prior to degradation and may be replaced by other, less demanding species.
In this work, we studied a mixed forest located in central Spain, in which seven tree species co-occur: Fagus sylvatica L., Quercus petraea (Matt.) Liebl., Quercus pyrenaica Willd., Ilex aquifolium L., Sorbus aucuparia L., Sorbus aria (L.) Crantz and Prunus avium (L.) L. The forest has been exploited since at least 1,460 for firewood and cattle grazing (Gil et al., 2009). This has resulted in an open woodland structure with scarce tree regeneration and scattered large, mature and multi-centennial trees regularly pollarded to obtain firewood. However, for conservation purposes, cattle grazing was banned in 1961 and the last logging operation was carried out in 1962 (López Santalla et al., 2003; Gil et al., 2009), favouring a rapid densification, increasing competition, and generating the current uneven-aged structure and species mixing. The forest was declared of national interest in 1974 and, together with other neighboring forests, it was declared UNESCO Biosphere Reserve since 2005. In this area, temperatures are increasing since the 1970s (Rubio-Cuadrado et al., 2020a), a trend that has been accentuated in the last decade (Rubio-Cuadrado et al., 2021a). In addition, the frequency of late-spring frosts is increasing (Sangüesa-Barreda et al., 2021). Taking into account the differences among species in shade and drought tolerance (Niinemets and Valladares, 2006), late frost susceptibility (Rubio-Cuadrado et al., 2021a), and competitive abilities (Hein and Dhôte, 2006; Manso et al., 2015), we aimed to study the effects of passive restoration on forest dynamics over the last decades. Our specific objectives were: (ⅰ) to characterize the current forest structure and composition, (ⅱ) to study the natural dynamics of the forest structure after its abandonment, (ⅲ) to identify the main structure drivers, and (ⅳ) to forecast a range of possible forest futures, including potential alterations to the current dynamics. Our hypotheses were:
H1 Due to the open forest structure at the time of the grazing and logging cessation, we expect that less shade-tolerant species (Q. petraea and especially Q. pyrenaica) initially increase their abundance until canopy closure is complete.
H2 After canopy closure, more shade-tolerant species (I. aquifolium and F. sylvatica) benefit the most, because (H2a) they are able to increase their recruitment and abundance in the juvenile classes, and (H2b) adult trees are better able to withstand increased competition with no loss of vigor and little mortality.
H3 Less drought-tolerant species, F. sylvatica and S. aucuparia, will tend to occupy higher altitudes and be more conditioned by soil depth due to ongoing climate change.
The study site, "El Hayedo de Montejo", is a mixed forest of 125 ha located at the "Sistema Central" mountain range in central Spain (41°07′ N, 3°30′ W; Fig. S1), between 1,250 and 1,500 m above sea level (a.s.l.). The site aspect is mainly east and the dominant species are Quercus pyrenaica, Fagus sylvatica and Quercus petraea. The medium-sized deciduous trees Sorbus aucuparia, Sorbus aria and Prunus avium, and the medium-sized evergreen tree Ilex aquifolium are present to a lesser extent, and other species such as Crataegus monogyna L., Taxus baccata L., Fraxinus excelsior L., Salix ssp. or Betula ssp. are marginally present. Except for the Mediterranean Q. pyrenaica, all other species are typical of temperate European climates and are in this forest near their southwestern distribution limit, especially in the case of F. sylvatica and Q. petraea (Fig. S2). Fagus sylvatica, Q. petraea, S. aucuparia, S. aria and P. avium, as well as most individuals of Q. pyrenaica originated from natural seedling in the studied forest, in contrast with the resprouting origin of Q. pyrenaica coppice stands so abundant in Central Spain (Valbuena-Carabaña et al., 2008). Ilex aquifolium originated mainly by stem sprouting, although individuals with seed origin also occur, so that vegetative reproduction is not a problem for genetic diversity in this forest (Vega et al., 2022).
All species show moderate to high shade tolerance. The most tolerant species is F. sylvatica, followed by I. aquifolium, P. avium, S. aria, S. aucuparia, Q. petraea and Q. pyrenaica. In terms of drought tolerance, these species show moderate to low tolerance. The most drought-tolerant species is Q. pyrenaica, followed by S. aria, I. aquifolium, Q. petraea, P. avium, F. sylvatica and S. aucuparia (Table S1) (Niinemets and Valladares, 2006; Rodríguez-Calcerrada et al., 2015).
The climate is Continental-Mediterranean, with 900 mm of mean annual rainfall and a prolonged dry period during July and August. The mean annual temperature is 9.5℃ (climate data from a meteorological station located in the forest and managed by the "Functioning of Forest Systems in a Changing Environment, FORESCENT" research group of the "Universidad Politécnica de Madrid"). The soil has been classified as humic cambisol and a horizon reaches 50 cm of depth on average, enabling considerable water storage during dry periods (Pardo et al., 1997).
We recovered the first inventory carried out in the study site, which dates back to 1909 when the forest was kept as a woodland pasture for livestock grazing (de Madariaga, 1909). This inventory consisted in a full-cover inventory in which only F. sylvatica trees were measured. More recently, three forest inventories were carried out in 1994, 2005 and 2015. In each inventory, 125 circular plots were systematically distributed on a square grid of 100 m × 100 m covering the entire forest area (Fig. S3). In the 1994 inventory, a plot radius of 12.6 m was used (plot area of 500 m2) for measuring the diameter at breast height (DBH) of trees with DBH greater than 5 cm. Fagus sylvatica, Q. petraea, Q. pyrenaica, I. aquifolium, S. aucuparia, S. aria and P. avium tree number per ha and basal area were assessed. In addition, to characterize the natural regeneration of the forest, a concentric subplot of 8 m radius was located within every plot to measure plant height and plant basal diameter of the saplings, ranging from newly germinated plants up to a DBH of 5 cm.
In the 2005 inventory, the plot radius was increased to 15 m (plot area of 707 m2), for measuring DBH of all trees with a minimum threshold of 5 cm, including snags and logs. Plots, as well as live trees, snags and stumps within the plots, were georeferenced. Vigor of all trees was noted as dead, for individuals with no living buds; very decaying, for individuals structurally damaged (i.e. decapitated and/or with dried main stem branches and/or severe rot infection) but with some living branches; decaying, for individuals with generalized top drying and/or some dead main branches and/or medium pathogen infection; and vigorous or very vigorous, for individuals with few or no signs of decay. Since 2005, plots were permanent, so the plots in the second and third inventories were at the same locations.
Finally, in the 2015 inventory, the plot design was based in two concentric circles of radius 10 and 15 m. All living and dead trees with a DBH between 5 and 20 cm were measured within the small circle and all living and dead trees with DBH larger than 20 cm were measured within the larger circle. As in the 2005 inventory, the vigor of all trees was recorded using the same scale. To characterize the natural regeneration of the forest, 40 plots among those 125 plots inventoried, were randomly selected and one belt transect (20 m × 4 m) was established along the North-South diameter of the plot. Species identification, plant height and plant basal diameter were recorded for each sapling (see Perea et al., 2020). Soil depth was also measured in the 125 plots. For this purpose, a 1-m long metal bar was hammered into four points in each plot. The bar was driven until there was a noticeable increase in penetration resistance, due to the bar reaching the lower horizon, which was highly compacted, with a high proportion of silts and clays, and often with the presence of disaggregated rock derived from the shale parent rock. For each plot, altitude, slope and aspect were obtained from a Digital Terrain Model (PNOA, Instituto Geográfico Nacional, Spain) with a spatial resolution of 25 m.
Additionally, two full-cover inventories in 2009 and 2019, evaluated all trees with DBH larger than 40 cm (estimated to be centennial). In the 2009 inventory, F. sylvatica, Q. petraea and Q. pyrenaica were evaluated, but in the 2019 inventory, only F. sylvatica were considered. The location, DBH and vigor status of the trees were measured. The vigor status was classified as dead (dead trees were only counted for F. sylvatica); very decaying; decaying; vigorous, for individuals with partial top drying and/or low pathogen affection; and very vigorous, for individuals with no sign of decay.
The sampling design was modified by using different plot sizes and spatial arrangement. These changes may affect the estimation error, however, due to the relatively high sampling intensity (one plot per hectare), the inventory result can be considered representative of the actual forest composition and therefore comparisons between inventories can be considered reliable, except for low abundance and heterogeneously distributed species and DBH classes.
To analyze the natural dynamics of the forest, values of tree number per ha (N) and basal area (BA) per diameter classes were calculated for each species in each inventory. To avoid DBH classes with too few trees, an increasing diameter class size was established. A total of 8 DBH classes were established: 5–10, 10–15, 15–20, 20–30, 30–40, 40–60, 60–80 and > 80 cm. In the case of regeneration, three height classes were established: < 0.3, 0.3–1.5 and > 1.5 m.
Differences between inventories were calculated using the Mann Whitney U test for independent samples and the Wilcoxon test for related samples (specifically to compare the 2005 and 2015 inventories, since both were measured at the same points). When the frequency of zeros was high in one of the groups (i.e., a diameter or height class is not represented by any tree in a high proportion of plots in one of the inventories), the differences between groups were analyzed with hurdle models using the pscl package (Zeileis et al., 2008) in R statistical software (R Development Core Team, 2022). The hurdle regression is a two-component model: a truncated count component is employed for positive counts, and a hurdle component models zero vs. larger counts. The year of inventory has been included as the only predictor in both components of the model. A negative binomial and a binomial distribution were employed for the count and hurdle components respectively. In those cases where the negative binomial distribution was not met by data, the count component of the hurdle model was replaced by the Mann Whitney U test.
Relationships between the different RE indexes obtained (see section 2.4 below), and between the mortalities of the different species, and the influence of stand and physiography on mortality and RE indexes, were analyzed by means of Spearman correlations. Data from all plots (one value per variable and plot) from the 2015 inventory are used for this purpose.
Hurdle regression models have been fitted in order to analyze the relationship between physiographic, soil and competition variables, as predictors, and BA and regeneration of each species and total (the aggregate of all the species), as response variables. As physiographic predictors in these models, altitude, slope and aspect of each plot (the latter variable transformed by means of the cosine to obtain values between +1 and −1, corresponding to north and south, respectively), were used. In addition, average soil depth per plot, and basal area of the different species per plot, as a proxy of competition, were used as predictors. In models of regeneration, the recruitment of species other than the one considered in the response variable was also included as predictors, to take into account the possible aggregation between species. Hurdle models were selected both to correctly address the high proportion of zeros in the response variables, and to be able to independently analyze the predictors that influence the two components of the models. Within the hurdle component, the relationships between the response variables and the BA and regeneration predictors indicate aggregation or rejection relationships among species, i.e., whether the increased presence of a species increases (positive coefficient) or decreases (negative coefficient) the probability of presence of the species considered in the response variable. The sign of any one of these relationships is independent of the competition effect in the count component, since, for example, an aggregation relationship may indicate competition for space. A negative binomial and a binomial distribution were employed for the count and hurdle components, respectively. When the response variable is aggregated (as occurs for total basal area or total regeneration, which result from the sum of the values of these variables for all the species), we adjusted negative binomial generalized linear models, using the MASS package in R (Venables and Riplley, 2002), because in these cases there were a low proportion of zero values. To identify the best-supported model, we constructed all possible combinations of alternative models from the full model (including all predictors), and then selected the best model by minimizing AIC. The existence of multicollinearity among explanatory variables was evaluated by calculating the variance inflation factor (VIF) using the packages performance (Lüdecke et al., 2021), for the hurdle models, and car (Fox and Weisberg, 2011) for the negative binomial generalized linear models. We used a threshold value of 4 for VIF to avoid the collinearity among predictors (Dormann et al., 2013). The goodness of fit of the models was assessed using the R2, understood as the square of the correlation between the actual values in the plots and those predicted by the model. For the hurdle models, R2 was calculated for the count component and the percentage of correct predictions was calculated for the hurdle component.
To evaluate the species dynamics, in particular the substitution of one species for a different one, we calculated the Species Replacement Indexes (REglobal; Ledo et al., 2009), which reflect the changes in species dominance based on the differences in diameter between pairs of trees belonging to the same or different species from a single inventory data. To obtain REglobal, the first step is to calculate the rglobal(Δ) function, which is based on the average conditional density (Stoyan and Stoyan, 1994). The rglobal(Δ) function averages a weighting function f(i, j), which assess the differences in composition, for all the pairs of trees i, j with a diameter difference Δ within the plot. In practice, f(i, j) is averaged over a bandwidth δ around Δ:
| (1) |
where n is the number of trees in a plot; Δ is the diameter difference; f(i, j) is the weighting function (that compares all pairs of trees in a plot for a given Δ), which takes value −1 when tree i (i = 1, 2, …, n) and tree j (j = i + 1) belong to the same species and 1 if the i species does not coincide with j; and the Iij (Δ ‒ δ, Δ + δ) function gives 1 if the diameter difference between tree i and j is within the (Δ ‒ δ, Δ + δ) interval or 0 if not (to include only those diameter differences that fall within the interval considered). The term Iij(Δ ‒ δ, Δ + δ) in the denominator ensures that rglobal(Δ) is unaffected by tree density within the plot. We used diameter difference intervals (Δ ‒ δ, Δ + δ) of 0.05 m width centered on Δ = 0.05, 0.10, 0.15, …, 0.80 m.
The r function can be calculated globally to assess any change in species composition (r(Δ)global) or, alternatively, it can be used to analyze the species A replacement with species B (r(Δ)A-B). In the case of r(Δ)A-B, f(i, j) takes value 1 if the tree with the smallest diameter is of the species B, and it takes value −1 when the smallest diameter is of the species A, and 0 if the pair of trees does not belong to the two considered species.
The value of the r(Δ)global function is close to −1 when, for a given Δ, we find conspecific pairs of individuals, and the value is close to 1 when trees belong to different species. In a similar way, the value of the r(Δ)A-B function is close to −1 when, for a given Δ, species A dominates among the trees with smaller diameters, and the value is close to 1 when species B dominates among the trees with smaller diameters (i.e. species B is replacing species A).
Finally, the Replacement Species Index (RE) is defined for each plot and inventory as the sum of the r(Δ) functions weighted by Δ:
| \mathrm{RE}=\sum\limits_{i=1}^t \Delta_i r\left(\Delta_i\right) | (2) |
where Δi is the center of the diameter difference interval i, i is the distinct diameter difference intervals present on the plot considered (normally only some of the preset intervals appear on each plot), and t is the largest diameter difference interval, which is set at 0.8 m, so that all those differences greater than 0.8 m were included in Δ = 0.8. This maximum was established to avoid giving excessive weight to larger trees, with low presence in the studied forest. We obtain the replacement index REglobal or REA-B, for r(Δ)global or r(Δ)A-B respectively. RE is not dimensionless and the result depends on the Δ units. In this study, Δ values used for weighting were measured in meters. Moreover, the absolute value of these indexes (REglobal and REA-B) increases with the width of range of diameter values in the plot.
Since the values across which these indexes move are not bounded, in order to interpret them reliably, it is necessary to analyze their evolution over time. An increasing value over time in REglobal indicates increasing different-species pairs of trees between smallest and largest diameter classes, whereas an increasing value over time in REA-B indicates that species B is becoming more abundant among the trees with smaller diameters (in the lower strata) in those areas where A was dominant. Note that REA-B = ‒ REB-A. To facilitate understanding of these indexes, some examples are provided in section 8.3 within the Supplementary Material.
According to the inventory of 2015, F. sylvatica and Q. pyrenaica were the most abundant species, with 30% and 27% of total trees (N) and 23% and 45% of the basal area (BA), respectively (Figs. 1 and S3, Tables S2 and S3). Fagus sylvatica was the most abundant species in the lower diameter classes and the second most abundant in the upper diameter classes, while Quercus pyrenaica was the most abundant species above 20 cm of DBH. Q. petraea and I. aquifolium showed a similar abundance, ca. 18%. However, whereas Q. petraea was homogeneously distributed across all diameter classes, I. aquifolium was mostly present in DBH classes below 30 cm. Finally, S. aucuparia, P. avium and S. aria were the least abundant species, with 3.2%, 2.8% and 0.4% of N respectively, and with trees only in the lower and intermediate diameter classes (< 40 cm). While P. avium and S. aria were homogeneously distributed among these diameter classes, most individuals of S. aucuparia belonged to the 5–10 and 10–15 cm diameter classes.
In the historical inventory carried out in 1909, the density of F. sylvatica was 49.2 trees·ha−1 and the BA 9.2 m2·ha−1 (Table S4). De Madariaga (1909) expressly noted that Q. petraea and Q. pyrenaica were almost all large, old, without regeneration, and considerably less abundant than F. sylvatica. In addition, I. aquifolium was abundant at that time.
In 1994, thirty years after grazing and logging cessation, Q. pyrenaica was the most abundant species, representing 28% of the total number of trees (Table S2). The diameter distribution showed in 1994 a gradual increase in the proportion of Q. pyrenaica from the highest diameter class to the 30–40 cm diameter class, where this species represented 29% and 68% of the trees, respectively (Fig. 2 and Table S2). Below this diameter class, the trend is reversed, with an abrupt fall in the proportion of Q. pyrenaica in the lowest diameter classes, where it represents 20% of the trees. In contrast, the proportions of F. sylvatica and especially Q. petraea showed the opposite distribution along the diameter classes.
Between 1994 and 2005, there was a generalized increase in density and BA (Fig. S4 and Table S5). However, from 2005 to 2015 the increase occurred only in the intermediate diameter classes, between 15 and 40 cm, with a decrease in the number of trees in the smallest diameter class (5–10 cm), which results in an increase in basal area between these two inventories, but a stagnation in density.
The long-term census data show a shift in species-specific relative abundance towards larger diameter classes between 1994 and 2015 due to the growth of trees between inventories (Fig. 2). The most interesting changes in relative densities between species occurred at the two extremes of diameter size classes. While in 1994 the proportion of F. sylvatica, Q. petraea and Q. pyrenaica was balanced for the smallest and largest diameter classes (32%, 18% and 20% of trees of 5–10 cm DBH vs 31%, 40% and 29 % of trees > 80 cm DBH belonged to these species, respectively; Table S2), in 2015 there was a significant increase in the proportion of F. sylvatica for the smallest diameter classes and the proportion of Q. pyrenaica for the largest diameter classes (Table 1). Indeed, while Q. petraea dominated the largest diameter class according to the 1994 inventory, when 40% of the trees belonged to this species, in 2015 this proportion dropped to 18%. In contrast, Q. pyrenaica increased from 29% to 55% in this diameter class, due to the incorporation of new trees from intermediate diameter classes. In the > 80 cm DBH class, the absolute density of Q. pyrenaica steadily increased by 112% between 1994 and 2015, while F. sylvatica (−4%) stagnated, and Q. petraea (−50%) decreased sharply (Table S5). It should be noted though, that these differences may partly be due to differences between the sampled locations in the 1994 and 2015 inventories, since it was not until the 2005 inventory that the permanent plots were installed. The different location of the plots may affect the proportion of trees in higher diameter classes due to their low density and uneven distribution in the forest. This seems to be the case for Q. petraea, where a drop in density in the largest diameter class was observed from 1994 to 2005 but did not change between 2005 and 2015 (Table S2). In the case of F. sylvatica, there was an increase in density in the largest diameter class between 1994 and 2005, possibly due to differences in the sampled area, but a sharp drop in density in 2015 due to mortality. Finally, the diameter class > 80 cm for Q. pyrenaica followed the same upward trend across the three inventories.
| Diameter class (cm) | F.sylvatica | Q.petraea | Q.pyrenaica | I.aquifolium | S.aucuparia | S.aria | P.avium | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 5–10 | 0.06+* | −0.03–* | 0.01 | 0.03 | –-0.06+*–x | 0.000 | −0.01 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 10–15 | 0.04+* | −0.08–* | 0.02 | −0.01 | 0.03+* | −0.002–* | 0.00 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 15–20 | 0.12+x | 0.04 | −0.12 | −0.07 | 0.00 | 0.001 | 0.03 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 20–30 | 0.09+x | 0.11+x | −0.15–*+x | −0.05–* | 0.00 | −0.005 | −0.01 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 30–40 | 0.05+x | 0.16+x | −0.15–*+x | −0.05–* | 0.00 | 0.007 | −0.01 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 40–60 | −0.04 | 0.06+x | 0.01+x | −0.01 | −0.01 | 0.000 | 0.00 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 60–80 | 0.00 | −0.11 | 0.12+x | 0.01 | 0.00 | 0.000 | −0.02 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| > 80 | −0.04 | −0.22 | 0.26+x | 0.00 | 0.00 | 0.000 | 0.00 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| Asterisks (*) indicate significant differences using Mann Whitney U test for the values larger than zero and symbols (X) indicate significant differences in the hurdle component of the hurdle regression (zero vs. larger values). Symbol (+) indicates greater presence (that is, a lower proportion of zeros) or greater relative density of the species considered (without considering the zeros), depending on the test considered. Symbol (−) indicates the opposite. When p-value is between 0.05 and 0.1 the asterisk/symbols has been underlined; otherwise p-value is less than 0.05. | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
DownLoad:
CSV
The relative density of big I. aquifolium trees decreased from 1994 to 2015 (Table 1). While in 1994 its relative density represented 19% and 12% of the trees in the diameter classes 20–30 and 30–40 cm, respectively; in 2015, it represented 14% and 8%, respectively, reflecting an increase in density over time lower than that of other species (Table S5). However, this species showed the second highest increase (after F. sylvatica) in relative density for the 5–10 cm DBH class. Sorbus aucuparia displayed a decrease in absolute (Table S5) and relative density (Table 1) in the smallest diameter class. In contrast, there was a significant increase in the relative density of S. aucuparia for the diameter class 10–15 cm, with the rest of the diameter classes remaining unchanged. Sorbus aria showed a decrease in relative density for diameter class 10–15 cm (Table 1), although in absolute terms density in S. aria increased (Table S5). Finally, P. avium showed no significant changes in relative density. Fagus sylvatica, I. aquifolium and S. aria were the only species that showed significant increases in (absolute) density in the smallest diameter class (5–10 cm; Table S5).
The replacement index (RE) values obtained from the 1994 inventory data showed a very similar species composition between the largest and smallest DBH. However, already in this year REQ. pyrenaica–F. sylvatica showed a higher proportion of F. sylvatica among the trees with small diameter (Table 2). In the 2005 inventory a clear differentiation of species according to the diameter class was observed, a trend that was more evident in the 2015 inventory. Fagus sylvatica was replacing Q. pyrenaica and, to a lesser extent, Q. petraea in stands where they were mixed. Quercus petraea was also replacing Q. pyrenaica where they grew together.
| Inventory | REglobal | REQ.pyrenaica – F. sylvatica | REQ.pyrenaica – Q. petraea | REF.sylvatica – Q. petraea |
| 1994 | −0.03 | 0.19 | −0.08 | 0.03 |
| 2005 | 0.21 | 0.65 | 0.13 | −0.23 |
| 2015 | 0.22 | 0.88 | 0.23 | −0.28 |
DownLoad:
CSV
In the 2015 inventory, four species showed a significant increase in density of recruits (compared to 1994; Fig. 3 and Table S6). Fagus sylvatica was the species with the highest density of regeneration (1,125 saplings·ha−1 considering the three types of regeneration together), and had increased its presence by 70% since 1994; Q. petraea, with 653 saplings per ha, increased by 150%; I. aquifolium, with 369 saplings per ha, increased 24%; and S. aria, with 38 saplings per ha, increased 64%. In contrast, Q. pyrenaica, with 459 saplings per ha, lost 52% saplings since 1994, when it was the species with most recruits; and S. aucuparia, with 66 saplings per ha, also decreased ca. 36%. Prunus avium had 347 saplings per ha, showing no variation since 1994. Considering only the largest height class, which is the most likely to survive the competition and growth up to the canopy, F. sylvatica represented 78% of the saplings, followed by Q. pyrenaica (12%), I. aquifolium (6%) and Q. petraea (4%) (Fig. 3d).
Fagus sylvatica showed a low mortality in DBH classes below 40 cm, and > 10% above this diameter class, reaching 37% for the largest diameter class (Table 3; note that the comparison was only between the 2005 and 2015 inventories, as the same permanent plots were measured in both years). On the contrary, Q. pyrenaica showed low mortality in diameter classes larger than 20 cm, but high mortality in the lower diameter classes, especially for the 5–10 cm diameter class. Quercus petraea is the species that showed, together with S. aria, the lowest mortality, being only important in the smallest diameter class. Ilex aquifolium and, especially, S. aucuparia and P. avium showed high levels of mortality in practically all diameter classes.
| Diameter class (cm) | F.sylvatica | Q.petraea | Q.pyrenaica | I.aquifolium | S.aucuparia | S.aria | P.avium |
| 5–10 | 6.62 (5.08%) | 6.88 (11.07%) | 20.12 (23.33%) | 6.62 (10.17%) | 13.24 (30.87%) | 0 (0%) | 2.04 (16.67%) |
| 10–15 | 1.53 (2.69%) | 1.53 (3.72%) | 6.11 (11.76%) | 5.09 (9.26%) | 1.27 (15.41%) | 0 (0%) | 2.55 (33.58%) |
| 15–20 | 1.02 (4.81%) | 1.02 (3.41%) | 3.82 (12.83%) | 4.07 (15.19%) | 0.76 (42.19%) | 0 (0%) | 0.51 (10.47%) |
| 20–30 | 0.34 (2.24%) | 0.11 (0.5%) | 2.49 (6.47%) | 1.92 (11.33%) | 0 (0%) | – | 0.11 (7.69%) |
| 30–40 | 0 (0%) | 0 (0%) | 0.23 (1.35%) | 0.45 (16.67%) | – | 0 (0%) | 0.11 (33.33%) |
| 40–60 | 0.57 (12.2%) | 0 (0%) | 0.34 (2.42%) | 0.11 (20%) | 0.11 (100%) | – | 0.11 (100%) |
| 60–80 | 0.23 (11.11%) | 0 (0%) | 0.11 (1.56%) | – | – | – | – |
| > 80 | 0.79 (36.84%) | 0 (0%) | 0 (0%) | – | – | – | – |
DownLoad:
CSV
All species, except for I. aquifolium and S. aria, showed a decrease in vigor in the lowest diameter classes (Fig. 4; note that this figure considers vigor and mortality in both 2005 and 2015). In the rest of the diameter classes, large differences were found among species. Fagus sylvatica showed a significant drop in vigor in upper diameter classes. While the number of very vigorous trees of F. sylvatica decreased by 9.5% in the full-cover inventory, and the number of vigorous trees decreased by 34.2%, the number of very decaying trees increased by 64.1% and mortality was close to 10%, reaching 14.4% for trees of F. sylvatica over 80 cm DBH (Table S7). On the contrary Q. petraea showed only slight (but significant; Fig. 4) declines in vigor below 20 cm DBH. Quercus pyrenaica showed the greatest loss of vigor, being generalized for all DBH classes, except for the highest (> 80 cm; Fig. 4). Ilex aquifolium showed considerable mortality in both 2005 and 2015 (Fig. 4), with significant increases in the intermediate diameter classes. Sorbus aucuparia and P. avium also showed declines in vigor in the intermediate diameter classes, although not significant due to the scarce presence of these species for these sizes.
Basal area of each species was inversely related to the BA of the rest of species, especially in the case of F. sylvatica, Q. petraea and Q. pyrenaica, as a result of competition for space (count component in Table 4). In the case of I. aquifolium and P. avium, the only species that competes with them for space was Q. pyrenaica (Table 4). Fagus sylvatica and Q. pyrenaica showed a rejection relationship with each other, while F. sylvatica showed a tendency to aggregate with Q. petraea (hurdle component in Table 4). Also I. aquifolium showed a tendency to aggregate with F. sylvatica and Q. petraea, while S. aucuparia showed rejection with F. sylvatica and Q. pyrenaica. Finally, BAother (result of the sum of the BA of I. aquifolium, S. aucuparia, S. aria, P. avium and other species with lesser presence in the forest) was directly related to the probability of occurrence of F. sylvatica, Q. petraea and P. avium (excluding P. avium from BAother when fitting P. avium model). The absence of predictors in the S. aria model cannot be interpreted as independence with respect to the variables considered; rather, it reflects the low abundance of the species.
| Basal area. Count componenta | ||||||||
| F.sylvatica | Q.petraea | Q.pyrenaica | I.aquifolium | S.aucuparia | S.aria | P.avium | All speciesb | |
| BAF.sylvaticac | −0.35 (0.16)* | −0.24 (0.08)** | ||||||
| BAQ.petraeac | −0.28 (0.14)* | −0.31 (0.13)* | ||||||
| BAQ.pyrenaicac | −0.79 (0.17)*** | −0.66 (0.18)*** | −0.65 (0.22)** | −1.16 (0.61)· | ||||
| Altitude | −0.58 (0.20)** | −0.36 (0.09)*** | −0.16 (0.06)** | |||||
| Slope | 1.03 (0.27)*** | 0.85 (0.42)* | ||||||
| Aspect | ||||||||
| Soil depth | 0.32 (0.14)* | 0.96 (0.43)* | 0.19 (0.06)** | |||||
| R2 | 0.42 | 0.14 | 0.37 | 0.47 | 0.31 | – | 0.17 | 0.11 |
| Basal area. Hurdle (zero) componenta | ||||||||
| F.sylvatica | Q.petraea | Q.pyrenaica | I.aquifolium | S.aucuparia | S.aria | P.avium | ||
| BAF.sylvaticac | −1.29 (0.32)*** | 0.39 (0.21)· | −1.38 (0.82)· | |||||
| BAQ.petraeac | 0.78 (0.33)* | 0.86 (0.25)*** | ||||||
| BAQ.pyrenaicac | −0.77 (0.27)** | −1.56 (0.74)* | ||||||
| BAotherc | 0.43 (0.26)· | 0.62 (0.27)* | 0.43 (0.22)· | |||||
| Altitude | −1.21 (0.3)*** | 0.92 (0.23)*** | −2.63 (0.49)*** | |||||
| Slope | −0.35 (0.21)· | 0.38 (0.21)· | ||||||
| Aspect | 0.84 (0.33)* | |||||||
| Soil depth | 0.37 (0.22)· | 0.58 (0.24)* | ||||||
| % correct | 79 | 71 | 83 | 73 | 86 | – | 67 | |
| a Predictors were standardized previously to fit the models. Note that positive coefficients in the hurdle component indicate that an increase in the factor increases the probability of a non-zero count. The goodness of fit of the model is shown by the R2 in the count component and the percentage of correct predictions in the hurdle component. b For the total basal area (considering all species) a negative binomial generalized linear model, instead of the hurdle model, was fitted due to the low proportion of plots where the response variable takes a value of 0. c BA is the basal area of the different species, with BAother being the total BA of the least abundant species (I.aquifolium, S.aucuparia, S.aria, P.avium and other species with lesser presence in the forest, but excluding, when apply, the species that correspond to the response variable, i.e., if the response variable is BAI.aquifolium, then BAother results from the sum of the basal areas of S.aucuparia, S.aria, P.avium and the other species with lesser presence in the forest). | ||||||||
DownLoad:
CSV
Fagus sylvatica and Q. pyrenaica were more abundant at lower elevations (Table 4) whereas Q. petraea mostly appeared at higher elevations. Slope had a positive direct relationship with the BA of I. aquifolium and S. aucuparia, but negative with F. sylvatica. Aspect was only related to Q. pyrenaica, a species that appears mostly in the northern orientation, contrary to what was expected since this is the most Mediterranean species in the forest. Finally, deeper soils were directly related to a greater presence of F. sylvatica, I. aquifolium and P. avium.
The occurrence and abundance of saplings taller than 1.5 m of the most shade tolerant species, F. sylvatica, I. aquifolium and P. avium, showed a positive relationship with the presence of the conspecific species, even if there is high competition (Table S8). Regeneration of the less shade-tolerant species did not occur or did not thrive under the conspecific canopy. In addition, there was an inverse relationship in the hurdle component between Q. pyrenaica regeneration and BA of F. sylvatica (i.e., the presence of F. sylvatica inhibits the occurrence of Q. pyrenaica regeneration). Altitude had a similar effect on the regeneration models as on the BA models, except for the case of F. sylvatica which, although present to a greater extent at lower elevations (Table 4), regenerated more at higher elevations (Table S8), indicating an expansion of beech towards higher altitudes. Regarding the other physiographic variables, slope was inversely related to the abundance of F. sylvatica regeneration, and the north-facing slope was positively related to the abundance of F. sylvatica and Q. petraea regeneration and to the presence of I. aquifolium. Counterintuitively, soil depth was inversely related to the abundance or presence of F. sylvatica, Q. petraea and S. aria regeneration (Table S8).
Mortality was related to increased competition, especially from the same species (Table S9). In areas with F. sylvatica, the high presence of species other than F. sylvatica reduces F. sylvatica mortality (Table S9). In general, mortality is lower in northerly aspects, in sites with deep soils, high altitudes and shallow slopes (Table S9).
In general, species replacement was favored in areas with higher competition and deep soils (see REglobal results in Table S10). Particularly, the replacement of Q. pyrenaica by F. sylvatica is favored by increased competition, whether of Q. pyrenaica or F. sylvatica (Table S10). However, high densities of Q. petraea and Q. pyrenaica in a given stand may delay their replacement by F. sylvatica and Q. petraea, respectively (Table S10).
At the end of the 1950s, the forest had an open structure over most of its area (Fig. S5), similar to the structure in 1909 (de Madariaga, 1909). After the prohibition of grazing, a densification of the forest took place (Figs. S3 and S4), from which Q. pyrenaica specially benefited. While in 1909 the most abundant species was F. sylvatica, in 1994 the most abundant species was Q. pyrenaica, with 28% of the trees and 42% of the basal area (Tables S2 and S3). The results of the full-cover inventory conducted in 2009 followed the same lines, and almost half of the trees were F. sylvatica (45%), even more so if we take into account the large number of dead trees of this species (see foot of Table S7). In this inventory, all trees of F. sylvatica, Q. petraea and Q. pyrenaica larger than 40 cm were sampled. Based on a previous study conducted at the same site and in which dominant trees were bored, it can be assumed that trees with DBH greater than 40 cm were generally established prior to 1961 (Rubio-Cuadrado et al., 2020a; see data used in Rubio-Cuadrado et al., 2020b). Therefore, we can consider this inventory of 2009 as an estimate of the abundance of these species prior to 1961.
The methodological differences between the inventories and the heterogeneous distribution of species in the forest (Fig. S3) make comparisons between inventories somewhat imprecise. However, this does not affect the interpretation of the trends found. Effectively, considering only the 1994 inventory, a gradual increase in the proportion of Q. pyrenaica is shown from the largest diameter class to the 30–40 cm diameter class, indicating an expansion of this species (Fig. 2). However, this trend is blurred in the 2015 inventory due to the mortality suffered by F. sylvatica in the higher diameter classes (Table 3 and Table S7).
There are numerous studies in the literature on Q. pyrenaica forests showing that, after being subjected to intensive management, once this is stopped, the forest remains with low vigor, low growth and lack of fruiting (Cañellas et al., 2004; Montes et al., 2004; Bravo et al., 2008; Salomón et al., 2017; Aldea et al., 2017; Moreno-Fernández et al., 2021). In all cases, however, they are coppice forests, a common structure in Q. pyrenaica forests, traditionally harvested for firewood and therefore with a strong root-to-shoot imbalance (Bravo et al., 2008; Salomón et al., 2013, 2016). In contrast, in the study area, Q. pyrenaica is a high forest of natural seed origin (Valbuena-Carabaña et al., 2008), with relatively high vigor (Fig. 4), growth (Rubio-Cuadrado et al., 2020a) and regeneration in 1994 (Fig. 3c).
Once grazing was banned, the open forest structure benefited Q. pyrenaica, the species with the lowest shade tolerance and higher light requirements of seedlings (Niinemets and Valladares, 2006; Rodríguez-Calcerrada et al., 2007, 2010). The light requirements of regeneration in this species leads to episodic expansion processes in open forest or in treeless areas, forming unimodal age cohorts, with tree recruitment decreasing drastically as the canopy closes (Camisón et al., 2015). In contrast, we found no evidence that Q. petraea, the second species present in the study area with the lowest shade tolerance, expanded after grazing ceased. In fact, its proportion remains more or less constant across all diameter classes (Fig. 2), which would indicate a constant and low level of recruitment, probably due to the higher water and nutrient requirements of Q. petraea saplings, in addition to differences in shade tolerance (Rodríguez-Calcerrada et al., 2008).
After canopy closure, the stand continued to densify until 2005, when the excess of competition leads to density stabilization (Fig. S4a) due to the balance between regeneration and mortality (Yaussy et al., 2013). At this point, the most competitive and shade-tolerant species, especially F. sylvatica (Manso et al., 2015; Table S1), out-persisted the rest of species. In temperate forests with full canopy, small forest gaps are filled with shade tolerant species that, either established and persisted as advance regeneration in the understory and regain growth as competition was released (Fig. 3; McCarthy, 2001), or colonize rapidly the small gaps created. In the studied forest, the RE index shows that during the first decades after the grazing ban, the same species present in the upper canopy regenerated under and around the adult trees (values close to zero for 1994 in Table 2). However, after stand densification, the more shade-tolerant species took advantage and regenerated under the canopy of the more intolerant species (see values for 2015 in Table 2), as shown by the increasing density in the lower diameter classes of F. sylvatica, I. aquifolium and S. aria, between 1994 and 2015 (Table S5). However, for I. aquifolium this increase was not reflected in the intermediate and upper diameter classes due to mortality (Table 3), probably caused by competition, although the influence of other factors, such as pests, could not be discarded (Castillo et al., 2009; McGregor et al., 2016). The relative density of F. sylvatica increased in all diameter classes below 40 cm (Table 1) and it is expected that in the future it will also increase in the higher diameter classes as the trees grow (Fig. 2). Increasing dominance of F. sylvatica would not be surprising in mixed beech-oak forests of central Europe (Petritan et al., 2014), but it was less expected in this rear-edge population, where summer drought is recurrent and occasionally severe. Fagus sylvatica is spreading in Central Spain despite being a species more sensitive to drought (Aranda et al., 2000, 2005; Piovesan et al., 2008) and late frosts (Rubio-Cuadrado et al., 2021b, 2021a) than oaks (Q. petraea and Q. pyrenaica) with which it competes, showing that the greater tolerance to low light of F. sylvatica is the determining factor explaining its dynamics in this forest (Rubio-Cuadrado et al., 2020a; de Tomás Marín et al., 2023). There are other factors that may also contribute to F. sylvatica expansion, such as the rapid growth in height of its saplings, which allows them to outgrow the saplings of the other species with which they compete and thus "win" the competition for light (Ligot et al., 2013; Van Couwenberghe et al., 2013), or the increasing populations of wild herbivores, preferring Quercus, Sorbus and Prunus over Fagus for browsing (Perea et al., 2020; Carpio et al., 2021).
The inventories showed a long-term decline of S. aucuparia, P. avium and Q. pyrenaica, and a stagnation of Q. petraea. The species with greater basal area is currently Q. pyrenaica; however, since 1994, it declined in relative density, especially in the lower-intermediate diameter classes (Fig. 2), pointing to lower establishment and loss of vigor (Figs. 3 and 4; see also Petritan et al., 2014). Sorbus aucuparia, P. avium and Q. petraea showed increases in density in the intermediate diameter classes, but decreases in the upper (in S. aucuparia and P. avium but not in Q. petraea; see differences between 2005 and 2015 in Table S2) and lower diameter classes (Table S5), with high mortality rates in S. aucuparia and P. avium (Table 3). Again, this is possibly related to increased competition, relatively low shade tolerance of these species (Table S1), and their progressive competitive exclusion by F. sylvatica and other trees dominating the canopy.
The most important factor determining the dynamics in the studied forest is the shade tolerance of the species, as shown by the expansion of the shade intolerant Q. pyrenaica when grazing was banned and the later expansion of the shade tolerant species, specially F. sylvatica, as canopy closed. Among the other secondary factors that influence the stand dynamic and structure, soil depth, slope and altitude stand out. Species distributions may be strongly driven by soil conditions (Diekmann et al., 2015; Walthert and Meier, 2017). In the studied forest, due to the relationship between soil depth, slope and water holding capacity and between climate and altitude, the less drought tolerant species, especially F. sylvatica, were more abundant and/or regenerated better in areas with gentle slopes, northerly aspects, deep soils and high altitude (Table 4, Tables S8 and S9; de Tomás Marín et al., 2023). Therefore, where these conditions are met, the replacement of less shade tolerant species with more shade tolerant species was accelerated (Table S10).
Ilex aquifolium tends to appear mixed with F. sylvatica and Q. petraea and, to a lesser extent, with P. avium (Table 4; note that I. aquifolium is, among the species included in the BAother variable, the most abundant, with the greatest BA, and therefore, with the largest weight in this variable). This is probably the result of facilitation between species, an effect already noted for I. aquifolium, whose seeds are dispersed by birds perching on trees and herbivores searching refuge under the canopy (Arrieta and Suárez, 2001, 2005).
Historical management is another important factor in explaining stand structure and dynamics. Fagus sylvatica is the dominant species with the highest mortality and loss of vigor in all the upper diameter classes (above 40 cm; Table 3 and Table S7). Indeed, de Madariaga (1909) counted 2,384 F. sylvatica trees over 40 cm DBH, most of which were > 70 cm DBH (Table S4), while in 2009 only 546 F. sylvatica trees over 40 DBH were counted. This situation could be explained by old age itself and, besides, traditional tree uses, as pruning, pollarding and topping for decades [potentially centuries (Pardo and Gil, 2005)] and competition from their own recruits. Nowadays, these centennial trees often show stem rots and symptoms of decay.
Regarding drought sensitivity, one may think that, being F. sylvatica the most sensitive species of the studied forest and being this forest at the meridional limit of the distribution of F. sylvatica, this species would grow under limiting climatic conditions and would suffer from low resilience and high mortality. This is what occurs in fir (Abies alba) forests and in some Scots pine populations located in their southernmost distribution limit (Oliva and Colinas, 2007; Giuggiola et al., 2010; Gazol et al., 2020). However, previous studies carried out in the same forest on trees established after the cessation of grazing, have shown higher vigor of F. sylvatica than other coexisting dominant species (Rubio-Cuadrado et al., 2020a), possibly due to drought adaptations specific to this F. sylvatica population (Aranda et al., 2015). This result is consistent with current expansion of F. sylvatica in other Spanish populations (Sánchez de Dios et al., 2016, 2020) and the high resilience found at the southern distribution limit of the species (Tegel et al., 2014; Cavin and Jump, 2017; Muffler et al., 2020). Collectively, these results do not corroborate the hypothesis that F. sylvatica grows under limiting conditions in rear edge (warm, drought-prone) populations.
The increasing dominance of F. sylvatica, which is the only species with effective regeneration (Fig. 3 and Table S6) and will progressively outcompete other species (Table 2), points to lower species diversity in the future. This advance regeneration will reach the canopy when competition release occurs, mostly when small canopy gaps open (McCarthy, 2001). This expansion of F. sylvatica at the expense of the rest of the species will likely occur in areas with gentle slope and deep soil, as well as at low altitudes, where this species is most abundant (Table 4), but also at higher elevations (Table S8), where F. sylvatica is expected to replace progressively the other species (mainly Q. pyrenaica; Table 2), either due to better climatic conditions or to the lesser competition of these species.
Quercus petraea is, among the dominant species, the only one with low mortality, except in the lower diameter class (Table 3) where it shows a loss of vigor and a decrease in density, both in absolute and relative terms (Table 1, Table S5 and Fig. 4). Its lower regeneration capacity under conditions of high competition (Fig. 3), and its higher browsing sensitivity (Perea et al., 2020), will lead to its decline in the long term, as has been described in other stands throughout Europe where this species is mixed with F. sylvatica (Petritan et al., 2017; Rubio-Cuadrado et al., 2018).
Under the Q. pyrenaica dominant stratum often appears F. sylvatica and, to a lesser extent, Q. petraea (Table S2), even in those sites with high basal areas of Q. pyrenaica (Table S10). In addition, Q. pyrenaica has suffered a loss of vigor in all diameter classes (Fig. 4). Although Q. pyrenaica is able to regenerate effectively in situations of moderate competition (Gea-Izquierdo et al., 2015), it has undergone a loss of weight in the lower diameter classes with respect to F. sylvatica and Q. petraea (Fig. 2), a decrease that will be transferred to the intermediate diameter classes in the short and medium term.
The tendency of I. aquifolium to appear mixed with F. sylvatica could lead to an expansion of I. aquifolium together with F. sylvatica (Table 4); however, I. aquifolium rarely reaches the intermediate and higher diameter classes due to its relative small size at maturity and its low growth rate (Rodríguez-Calcerrada et al., 2019), which, along with its high mortality (Table 3), led to a relative loss of weight in these diameter classes (Table 1). Apart from that, I. aquifolium can form very dense, monospecific patches of undergrowth, as in fact occurs in some parts of the study area (Vega et al., 2022), where high competition for light prevents other species from regenerating (Arrieta and Suárez, 2001).
Sorbus aucuparia and P. avium, similar to Q. pyrenaica, will gradually lose importance due to their lower recruitment (Fig. 3) and high mortality (Table 3). In the long term, Q. pyrenaica, S. aucuparia, P. avium and, to a lesser extent, Q. petraea, will be confined to areas with shallow soils and steeper slopes. The heterogeneous distribution of species throughout the forest (Fig. S3) and the high longevity of the dominant species will prevent F. sylvatica from dominating the entire forest in the medium term. However, in the absence of disturbances or drastic climatic changes, F. sylvatica will eventually dominate in the long term. The higher sensitivity of F. sylvatica to drought and late frosts, together with the increasing frequency of these climatic events, may lead to a loss of resilience of the forest, as already observed in the higher diameter classes (Table 3 and Table S7), more affected by senescence and anthropic uses. Analyses of forest dynamics in areas with high past anthropogenic pressure are essential, considering that forest resilience is not necessarily linked to the natural dynamics of the ecosystem (Buma and Wessman, 2013). In our case, if the aim is to maintain the current mix of species in the forest, forest management will have to favor the regeneration of the species that are currently declining and facilitate the growth of these species in the understory. On the other hand, managers will have to reduce competition on the singular, centennial trees that may see their vigor reduced and die if the competition continues to increase. Anyway, further studies are needed to study the causes behind the mortality of F. sylvatica trees in the higher diameter classes.
The study site is a mixed forest formed mostly by species typical of temperate European climates that are found here near their lower latitude (rear edge) end of the species distribution. However, the location of the forest in an area of Mediterranean influence does not seem to determine its dynamics, with drought resistance being apparently a minor factor influencing the developmental trajectory of the forest. On the contrary, it is the successional character of the species, their relative capacity to endure low-light conditions and compete when overtopped by mature trees, that seems to strongly determine the dynamics of the forest. The results showed that the abandonment of grazing and cutting benefitted the expansion of Q. pyrenaica, thanks to the lower light sensitivity of their recruits. However, once the stand densified, more shade-tolerant species, especially F. sylvatica, outcompeted the others and regenerated more successfully. If this pattern of forest dynamics continues, F. sylvatica will gradually overcome the rest of species, which will be confined to areas with shallow soils and steeper slopes. Mature, multi-centennial individuals of F. sylvatica will likely disappear soon due to the special conditions of low density and intense management in which they grew.
Álvaro Rubio-Cuadrado: Writing – review & editing, Writing – original draft, Methodology, Investigation, Formal analysis, Conceptualization. Iciar Alberdi: Writing – review & editing, Funding acquisition. Isabel Cañellas: Writing – review & editing, Funding acquisition. Fernando Montes: Writing – review & editing, Software, Methodology, Formal analysis. Jesús Rodríguez-Calcerrada: Writing – review & editing, Funding acquisition, Conceptualization. Rosana López: Writing – review & editing, Funding acquisition. Guillermo G. Gordaliza: Investigation, Data curation. María Valbuena-Carabaña: Investigation, Funding acquisition, Data curation. Nikos Nanos: Investigation. Ramón Perea: Writing – review & editing, Investigation, Funding acquisition, Data curation. Luis Gil: Writing – review & editing, Supervision, Resources, Funding acquisition, Conceptualization.
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Iciar Alberdi and Isabel Cañellas report financial support provided by Horizon Europe. María Valbuena Carabaña and Luis Gil report financial support provided by Community of Madrid, Spain. Ramon Perea reports financial support provided by Community of Madrid, Spain. Rosana Lopez, Jesús Rodríguez Calcerrada and Álvaro Rubio Cuadrado report financial support provided by Spanish Ministry of Science and Innovation, Spain. Ramon Perea reports financial support was provided by Spanish Ministry of Science and Innovation, Spain.
IA and IC acknowledges support by project SUPERB H2020 (Systemic solutions for upscaling of urgent ecosystem restoration for forest related biodiversity and ecosystem services). MVC and LG acknowledges support by project P2013/MAE-2760 (Autonomous Community of Madrid). RL, JRC, ARC, acknowledges support by project PID2019-107256RB-I00 (Spanish Ministry of Science and Innovation) and project FAGUS by the Comunidad de Madrid through the call Research Grants for Young Investigators from Universidad Politécnica de Madrid. RP acknowledges support by projects 9OHUU0-10-3L226X (Autonomous Community of Madrid), RTI2018-094202-BC21 and RTI2018-094202-A-C22 (Spanish Ministry of Science and Innovation).
Supplementary data to this article can be found online at https://doi.org/10.1016/j.fecs.2024.100197.
Figures(4) / Tables(4)
Supported by: Beijing Renhe Information Technology Co. Ltd. Technical support: info@rhhz.net
| Diameter class (cm) | F.sylvatica | Q.petraea | Q.pyrenaica | I.aquifolium | S.aucuparia | S.aria | P.avium | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 5–10 | 0.06+* | −0.03–* | 0.01 | 0.03 | –-0.06+*–x | 0.000 | −0.01 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 10–15 | 0.04+* | −0.08–* | 0.02 | −0.01 | 0.03+* | −0.002–* | 0.00 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 15–20 | 0.12+x | 0.04 | −0.12 | −0.07 | 0.00 | 0.001 | 0.03 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 20–30 | 0.09+x | 0.11+x | −0.15–*+x | −0.05–* | 0.00 | −0.005 | −0.01 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 30–40 | 0.05+x | 0.16+x | −0.15–*+x | −0.05–* | 0.00 | 0.007 | −0.01 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 40–60 | −0.04 | 0.06+x | 0.01+x | −0.01 | −0.01 | 0.000 | 0.00 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 60–80 | 0.00 | −0.11 | 0.12+x | 0.01 | 0.00 | 0.000 | −0.02 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| > 80 | −0.04 | −0.22 | 0.26+x | 0.00 | 0.00 | 0.000 | 0.00 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| Asterisks (*) indicate significant differences using Mann Whitney U test for the values larger than zero and symbols (X) indicate significant differences in the hurdle component of the hurdle regression (zero vs. larger values). Symbol (+) indicates greater presence (that is, a lower proportion of zeros) or greater relative density of the species considered (without considering the zeros), depending on the test considered. Symbol (−) indicates the opposite. When p-value is between 0.05 and 0.1 the asterisk/symbols has been underlined; otherwise p-value is less than 0.05. | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
DownLoad:
CSV
| Inventory | REglobal | REQ.pyrenaica – F. sylvatica | REQ.pyrenaica – Q. petraea | REF.sylvatica – Q. petraea |
| 1994 | −0.03 | 0.19 | −0.08 | 0.03 |
| 2005 | 0.21 | 0.65 | 0.13 | −0.23 |
| 2015 | 0.22 | 0.88 | 0.23 | −0.28 |
DownLoad:
CSV
| Diameter class (cm) | F.sylvatica | Q.petraea | Q.pyrenaica | I.aquifolium | S.aucuparia | S.aria | P.avium |
| 5–10 | 6.62 (5.08%) | 6.88 (11.07%) | 20.12 (23.33%) | 6.62 (10.17%) | 13.24 (30.87%) | 0 (0%) | 2.04 (16.67%) |
| 10–15 | 1.53 (2.69%) | 1.53 (3.72%) | 6.11 (11.76%) | 5.09 (9.26%) | 1.27 (15.41%) | 0 (0%) | 2.55 (33.58%) |
| 15–20 | 1.02 (4.81%) | 1.02 (3.41%) | 3.82 (12.83%) | 4.07 (15.19%) | 0.76 (42.19%) | 0 (0%) | 0.51 (10.47%) |
| 20–30 | 0.34 (2.24%) | 0.11 (0.5%) | 2.49 (6.47%) | 1.92 (11.33%) | 0 (0%) | – | 0.11 (7.69%) |
| 30–40 | 0 (0%) | 0 (0%) | 0.23 (1.35%) | 0.45 (16.67%) | – | 0 (0%) | 0.11 (33.33%) |
| 40–60 | 0.57 (12.2%) | 0 (0%) | 0.34 (2.42%) | 0.11 (20%) | 0.11 (100%) | – | 0.11 (100%) |
| 60–80 | 0.23 (11.11%) | 0 (0%) | 0.11 (1.56%) | – | – | – | – |
| > 80 | 0.79 (36.84%) | 0 (0%) | 0 (0%) | – | – | – | – |
DownLoad:
CSV
| Basal area. Count componenta | ||||||||
| F.sylvatica | Q.petraea | Q.pyrenaica | I.aquifolium | S.aucuparia | S.aria | P.avium | All speciesb | |
| BAF.sylvaticac | −0.35 (0.16)* | −0.24 (0.08)** | ||||||
| BAQ.petraeac | −0.28 (0.14)* | −0.31 (0.13)* | ||||||
| BAQ.pyrenaicac | −0.79 (0.17)*** | −0.66 (0.18)*** | −0.65 (0.22)** | −1.16 (0.61)· | ||||
| Altitude | −0.58 (0.20)** | −0.36 (0.09)*** | −0.16 (0.06)** | |||||
| Slope | 1.03 (0.27)*** | 0.85 (0.42)* | ||||||
| Aspect | ||||||||
| Soil depth | 0.32 (0.14)* | 0.96 (0.43)* | 0.19 (0.06)** | |||||
| R2 | 0.42 | 0.14 | 0.37 | 0.47 | 0.31 | – | 0.17 | 0.11 |
| Basal area. Hurdle (zero) componenta | ||||||||
| F.sylvatica | Q.petraea | Q.pyrenaica | I.aquifolium | S.aucuparia | S.aria | P.avium | ||
| BAF.sylvaticac | −1.29 (0.32)*** | 0.39 (0.21)· | −1.38 (0.82)· | |||||
| BAQ.petraeac | 0.78 (0.33)* | 0.86 (0.25)*** | ||||||
| BAQ.pyrenaicac | −0.77 (0.27)** | −1.56 (0.74)* | ||||||
| BAotherc | 0.43 (0.26)· | 0.62 (0.27)* | 0.43 (0.22)· | |||||
| Altitude | −1.21 (0.3)*** | 0.92 (0.23)*** | −2.63 (0.49)*** | |||||
| Slope | −0.35 (0.21)· | 0.38 (0.21)· | ||||||
| Aspect | 0.84 (0.33)* | |||||||
| Soil depth | 0.37 (0.22)· | 0.58 (0.24)* | ||||||
| % correct | 79 | 71 | 83 | 73 | 86 | – | 67 | |
| a Predictors were standardized previously to fit the models. Note that positive coefficients in the hurdle component indicate that an increase in the factor increases the probability of a non-zero count. The goodness of fit of the model is shown by the R2 in the count component and the percentage of correct predictions in the hurdle component. b For the total basal area (considering all species) a negative binomial generalized linear model, instead of the hurdle model, was fitted due to the low proportion of plots where the response variable takes a value of 0. c BA is the basal area of the different species, with BAother being the total BA of the least abundant species (I.aquifolium, S.aucuparia, S.aria, P.avium and other species with lesser presence in the forest, but excluding, when apply, the species that correspond to the response variable, i.e., if the response variable is BAI.aquifolium, then BAother results from the sum of the basal areas of S.aucuparia, S.aria, P.avium and the other species with lesser presence in the forest). | ||||||||
DownLoad:
CSV
Email Alerts
RSS Feeds