Competition overrides climate as trigger of growth decline in a mixed Fagaceae Mediterranean rear-edge forest

Abstract

Key message

In recent decades, there has been a decline in growth in a rear-edge broadleaf forest of Fagus sylvatica , Quercus petraea , and Q uercus pyrenaica . Although temperatures have been rising due to climate change, the observed decline in growth was mainly attributed to increased density and competition between trees since the cessation of traditional uses such as logging in the 1960s.

Context

In recent decades, two major factors have influenced tree growth in many forests: climate warming, which is associated with aridification and negative growth trends in many Mediterranean forests, and abandonment of forest management, resulting from forest policy in conjunction with rural depopulation in Europe, often leading to an increase in competition and a decrease in growth.

Aims

Here, we study the growth trends in a mixed forest of Fagus sylvatica, Quercus petraea, and Quercus pyrenaica, where the abandonment of traditional uses in the 1960s has been followed by an increase in tree density. In this forest, both F. sylvatica and Q. petraea reach their south-westernmost limits of distribution.

Methods

Using dendrochronological methods and growth modeling, we assess the importance of climate warming on the shifts in competitive growth advantage of these three coexisting tree species and the relative importance of climate and competition on growth trends.

Results

Q. petraea and especially F. sylvatica showed a favorable evolution of their competitive capacity, despite the increase in temperatures that has occurred in the area in recent decades. F. sylvatica presented the lowest sensitivity to climate.

Conclusion

Under the current climate and forest structure conditions, competition is the most limiting factor on tree growth for the two oak species.

1 Introduction

The last 30-year period has been the warmest on the Earth’s surface since 1850 (IPCC 2014). The Mediterranean region has been severely affected by this warming, which has brought about aridification of extensive drought-prone areas (Jacob et al. 2014). Climate change will particularly affect temperate and Mediterranean forests located at the equatorial limit (rear edge) of the tree species distribution (Hampe and Petit 2005). This is the case of certain hardwood tree species (oaks, beech) which form mixed forests in mountain areas of Spain under a typical sub-Mediterranean influence (Rubio-Cuadrado et al. 2018c), characterized by warmer and drier conditions in summer (Giorgi and Lionello 2008).

Climate change is already affecting forest growth and productivity (Reyer et al. 2014). Negative growth trends associated with climate warming have previously been described in several Mediterranean conifer forests (Macias et al. 2006; Sarris et al. 2007, 2011; Gea-Izquierdo et al. 2014; Camarero et al. 2015), but temperate oak and beech forests with recurrent summer droughts typical of the Mediterranean climate are also negatively impacted (Piovesan et al. 2008; Colangelo et al. 2017; Dorado-Liñán et al. 2017a). These negative growth trends are forecasted to increase over the course of the twenty-first century (Sánchez-Salguero et al. 2017). Indeed, the Iberian Peninsula has been described as the most sensitive region, and that which is most likely to suffer the largest decreases in forest productivity throughout the continent (Morales et al. 2007; Reyer et al. 2014).

It is unclear how mixed forests subjected to summer drought will be able to cope with the negative effects on growth of warmer and drier climate conditions. On the one hand, species mixing may improve the capacity of the forest to recover after disturbances (i.e., larger resilience; Pretzsch et al. 2013). This is due to resource use efficiency, which should enhance tree growth and stand productivity as a result of niche complementarity and facilitation between coexisting tree species (Richards et al. 2010). On the other hand, climate warming may favor species which are more adapted to drier and warmer conditions over typical temperate species (Rubio-Cuadrado et al. 2018b).

In addition to climate change, land use changes have also affected growth trends in European forests (Linares et al. 2009; Camarero et al. 2011). In recent decades, forest policy in Europe has moved towards minimal intervention or non-management for several different reasons: aiming to increase adaptability and resilience (Brang et al. 2014), for conservation purposes or due to the low economic value of forest resources (Moreno-Fernández et al. 2016). However, most European forests have a long history of anthropogenic transformations (Valbuena-Carabaña et al. 2010; Schoolman et al. 2018) and are composed of relatively young stands. In addition, during recent decades, forests have expanded into new areas, especially in Europe, mainly due to land use changes and afforestations, as well as to climate warming, all these new woodlands also being at early successional stages (Álvarez-Martínez et al. 2014; Penniston and Lundberg 2014; Price et al. 2017). Climate and land use changes have increased the frequency of disturbances (e.g., forest fires), thus preventing forests from reaching the latter serial stages (McLauchlan et al. 2020). In all these stands, natural dynamics has often led to an increase in tree density and tree-to-tree competition over recent decades, which can result in decreased tree growth (Voelker et al. 2008; Linares et al. 2010; Ruiz-Benito et al. 2013; Hosseini et al. 2018).

In this study, we analyze the growth trends in a mixed beech-oak forest (“El Hayedo de Montejo”; Fig. 1) dominated by Fagus sylvatica L., Quercus petraea (Matt.) Liebl., and Quercus pyrenaica Willd., with a long history of anthropogenic intervention (Pardo and Gil 2005). The studied forest harbors one of the southernmost (rear edge) and westernmost populations of F. sylvatica and Q. petraea in Europe, while in the case of Q. pyrenaica, it is the core of its distribution range (Pardo et al. 2004). Since the 1960s, traditional uses (logging, grazing) have ceased in the study area, and this has promoted forest recruitment and densification, evolving from open woodland to a dense forest (Gil et al. 2009). Previous studies conducted at this site have shown a general decline in growth of the dominant trees, attributed to temperature increase over recent decades (Dorado-Liñán et al. 2017a). The effect of competition on dominant trees is often neglected in dendroecological analyses, but substantial increases in density may reduce the resources available not only to dominated or suppressed trees but also to dominant trees, leading to growth decline. However, the complexity involved in disentangling the combined effects of climate change and the recent evolution of the forest, in particular the increase in competition, has meant that this situation has not yet been analyzed. The aim of this study is to evaluate the relative importance of climate and competition on the radial growth of dominant and co-dominant trees. Our specific objectives are as follows: (1) to study the evolution of the competitive growth advantage of the three coexisting species and (2) to analyze the relative importance of climate and competition on the radial growth of F. sylvatica, Q. petraea, and Q. pyrenaica in this mixed forest. We hypothesize that (1) the increase in temperature and frequency of drought events have a positive effect on the growth of the Mediterranean species (Q. pyrenaica) but have a negative impact on that of the temperate species (F. sylvatica and Q. petraea), (2) increase in competition has been the main limiting factor on the growth of the three species over recent decades, and (3) competition has a greater impact on Q. pyrenaica and Q. petraea than on the more competitive F. sylvatica. To fulfill the objectives and verify the hypothesis, we analyzed the relationships between climate, the tree-ring widths of the three species, and stand basal area. In order to evaluate the first hypothesis, we used climate-growth correlations and pointer year analysis. To evaluate the first and third hypotheses, the differences among the species, in terms of the long-term growth response to the evolution of site conditions, were analyzed using the competitive advantage index based on radial growth data (Rubio-Cuadrado et al. 2018b). Finally, to evaluate the second and third hypotheses, we analyzed the relative importance of competition and climate on basal area increment for the three species using mixed models.

Fig. 1

Map of the Iberian Peninsula with the distribution of Fagus sylvatica (blue color), Quercus petraea (red), and Quercus pyrenaica (green). Areas where at least two of the three species are mixed are represented in yellow. Sources: species distribution cartography has been elaborated by INIA-CIFOR from Third National Forest Inventory (https://www.miteco.gob.es/es/biodiversidad/servicios/banco-datos-naturaleza/informacion-disponible/ifn3.aspx), and Digital Terrain Model has been elaborated by © Instituto Geográfico Nacional

2 Materials and methods

2.1 Study area

The study site, “El Hayedo de Montejo,” is a mixed forest of 125 ha located in the Sistema Central mountain range, Central Spain (41° 07′ N; 3° 30′ W; Fig. 1), between 1250 and 1500 m a.s.l. The dominant species are Fagus sylvatica L. (hereafter beech), Quercus petraea (Matt.) Liebl. (hereafter oak), and the Mediterranean Quercus pyrenaica Willd., with a relative basal area of 27%, 23%, and 50% and a mean diameter at breast height of 16.5, 19.9, and 24.1 cm, respectively. This forest is close to the southernmost European distribution limit of the two temperate species (beech and oak). The beech and Q. petraea trees in the studied forest come from natural seeding, as do most of the Q. pyrenaica individuals, in contrast to coppice stands, which characterize Q. pyrenaica in Central Spain (Valbuena-Carabaña et al. 2008). The forest has historically been exploited for firewood and cattle grazing, which originated an open woodland structure of dispersed large, mature trees with scarce tree recruitment. However, for conservation purposes, cattle grazing has been forbidden in the forest since 1961, and the last logging operation was carried out in 1962 (López Santalla et al. 2003; Gil et al. 2009), favoring the increase in tree density and competition. In this regard, the inventories carried out in the forest in 1994, 2005, and 2015 (see Section 2.3 “Tree sampling and competition data”) gave an average basal area for these years of 20.9, 26.5, and 28.8 m² ha⁻¹, respectively, and an average density of 616, 952, and 910 trees ha⁻¹, respectively.

The site has a Mediterranean continental climate with a 900-mm mean annual rainfall and a prolonged dry period during July and August (Gil et al. 2009). The mean annual temperature is around 9.5 °C (climate data from the meteorological station of the “Natural Systems and Forest History” research group of the Universidad Politécnica de Madrid, located within the studied forest). The soil is classified as humic cambisol (Pardo et al. 1997), and the A horizon reaches an average depth of 50 cm which allows water to be stored during dry periods.

2.2 Climate data

Local climate data from the “Pantano El Vado” weather station (41° 0′ 13″ N and 3° 18′ 7″ W, 910 m a.s.l., AEMET, Spanish Meteorological Agency) located at ca. 21 km from the forest was used in this study. Missing monthly temperature and precipitation data as well as daily temperature data were estimated using linear regressions, which relate climatic variables from this station with the temperatures from “Puerto de Navacerrada” weather station (AEMET, 40° 47′ 35″ N and 4° 0′ 38″ W, 1894 m a.s.l.) and with precipitation data from “Arbancón” weather station (AEMET, 40° 58′ 0″ N and 3° 6′ 57″ W, 902 m a.s.l.) as follows:

$$ {\mathrm{MT}}_{\mathrm{Pantano}\ \mathrm{El}\ \mathrm{Vado}}=0.9786\ast {\mathrm{MT}}_{\mathrm{Puerto}\ \mathrm{de}\ \mathrm{Navacerrada}}+5.8616 $$

(1)

$$ {\mathrm{DT}}_{\mathrm{Pantano}\ \mathrm{El}\ \mathrm{Vado}}=0.847\ast {\mathrm{DT}}_{\mathrm{Puerto}\ \mathrm{de}\ \mathrm{Navacerrada}}+6.7033 $$

(2)

$$ {\mathrm{MP}}_{\mathrm{Pantano}\ \mathrm{El}\ \mathrm{Vado}}=1.2395\ast {\mathrm{MP}}_{\mathrm{Arbancon}}+3.9396 $$

(3)

where MT, MP, and DT and are the monthly data for temperature and precipitation and the daily data for temperatures, respectively. R² = 0.96 for eq. 1, R² = 0.83 for eq. 2 and R² = 0.76 for eq. 3. For the missing daily precipitation data from the “Pantano El Vado” local station, we used the daily data from “Arbancón” weather station (directly, without transforming) due to the low R² (0.13) regression line that relates both data series given the high spatial variability intrinsic to this variable.

In addition, to assess drought intensity, we used the Standardised Precipitation-Evapotranspiration Index (SPEI) (Vicente-Serrano et al. 2010a, b) calculated from the same local data using the Thornthwaite equation to estimate potential evapotranspiration within the SPEI package in R (Beguería and Vicente-Serrano 2017). The SPEI drought index considers the effect of temperature on evapotranspiration rates and the cumulative water deficit. High and low SPEI values correspond to wet and dry conditions, respectively. The period from the arrival of water input (precipitation) to availability of a given usable resource (water absorption by the roots) differs considerably. Thus, the time scale over which water deficits accumulate becomes extremely important to calculate the SPEI. For this reason, we used a range of time scales (from 1 to 20 months) during which the water deficit and surplus are accumulated (Vicente-Serrano et al. 2010b).

2.3 Tree sampling and competition data

Three forest inventories were carried out in 1994, 2005, and 2015 in the studied forest. In each inventory, 125 circular plots systematically distributed on a square grid of 100 m × 100 m were established covering the entire forest area. Permanent plots have been established since 2005, so the plots in the second and third inventories were at the same locations, whereas in 1994 the locations were different, although the area covered and the sampling intensity were the same. Based on the 2015 inventory data, all plots with mixed forest (presence of at least two of the three studied species) were selected for the extraction of two cores at 1.3 m (always perpendicular to the maximum slope) from a dominant or co-dominant tree of each species located nearest to the plot center using a Pressler increment borer. The diameter at breast height (DBH) of each tree bored was measured at 1.3 m with a caliper. Sampling took place during autumn of 2018. At least 20 trees per species, distributed across 32 plots, were sampled.

The total stand basal area (BA) and the fraction of it belonging to each of the species were calculated in 12.6-m radius plots for the 1994, 2005, and 2015 inventories. The area of the plot with this radius is slightly larger than the crown projection of these species, in accordance with the DBH-crown diameter relationship of the Second Spanish National Forest Inventory (DGCONA 1998), so the BA estimated can be considered an indicator of stand level competition for resources. Spatio-temporal Universal Kriging (UK) was used to interpolate the total BA and spatio-temporal Ordinary Kriging (OK) to interpolate the proportion of each species for the 1994–2015 period from the BA measured in the inventory plots. As vegetation indices can indicate changes in stand density (Gómez et al. 2012), normalized difference vegetation index (NDVI) averaged from the available Landsat images for July–August of each year was used as an auxiliary variable in the UK interpolation of BA. These models interpolate the BA value and the proportion of each species at the plot locations for each year of the 1994–2015 period using the spatial and temporal autocorrelation of the variables to determine the optimal kriging weights. The low annual growth rates of the three analyzed species in the study area along with the species regeneration and development processes result in a high spatio-temporal autocorrelation of BA at the distance and time lags between measurements (100 m distant sampling plots and 10 years between inventories), providing consistent interpolated values for the 1994–2015 period. The spatio-temporal autocorrelation was modeled through the generalized product-sum variogram model (De Iaco et al. 2002). The OK variogram parameters were estimated using the weighted least squares (Cressie 1985), whereas the UK variogram parameters and mean function coefficients were estimated through the iteratively reweighted generalized least squares method proposed by Neuman and Jacobson (1984). Cross-validation was carried out to verify the unbiasedness of the estimated values and the ratio between the estimated kriging variance and the mean of the cross-validation squared residuals. Competition data are available in Rubio-Cuadrado et al. (2020).

2.4 Radial growth data

Cores of the sampled trees were mounted on wooden supports and carefully sanded until tree rings were clearly visible. After visual cross-dating, tree-ring widths (TRW) were measured to the nearest 0.01 mm using the semi-automatic LINTAB measuring device with the TSAP-Win software (RINNTECH, Heidelberg, Germany). Cross-dating was then further verified using the COFECHA program (Holmes 1997). A beech sample, with a very different age (129 years) and growth trend and very low correlation with the mean species series (close to 0), was discarded from further analysis. TRW data are available in Rubio-Cuadrado et al. (2020).

We calculated the mean growth series for each individual, and these series were subsequently transformed into basal area increments (BAI) as this variable is better than tree-ring width for capturing growth trends and accounting for the increase in size and age of stems (Biondi and Qeadan 2008). The BAI was calculated as follows:

$$ BAI=\pi \left({r}_t^2-{r}_{t-1}^2\right) $$

(4)

where r_t and r_t-1 are the stem radii at the end and the beginning of a given annual ring increment corresponding to rings formed in t and t-1 years, respectively.

When the cores did not contain the pith, the length of the missing radius was estimated to calculate the BAI as the distance to the geometric center of the tree, assuming concentric growth. Estimation was carried out by subtracting the corresponding cored length to the measured radius in the field, previously subtracting an estimated width of the bark. Bark thickness was estimated based on tree diameter using equations obtained from the data of the Second Spanish National Forest Inventory (DGCONA 1998). The functions used to estimate the bark thickness of F. sylvatica (C_Fs), Q. petraea (C_Qpe), and Q. pyrenaica (C_Qpy) are the following:

$$ {C}_{Fs}=\kern0.5em -0.00001\ {\mathrm{DBH}}^2+0.0284\ \mathrm{DBH}+0.3545 $$

(5)

$$ {C}_{Qpe}=1.4407\ {\mathrm{DBH}}^{0.4071} $$

(6)

$$ {C}_{Qpy}=7.4422\ \ln\ \left(\mathrm{DBH}\right)\hbox{--} 23.845 $$

(7)

The final species chronologies were built by averaging all individual BAI measurements of the same species. The statistical quality of each chronology was checked via expressed population signal (EPS) (Table 1). The sample depth is considered adequate when EPS > 0.85 (Speer 2012).

Table 1 Description of dendrochronological variables of cored trees

Pointer years (years with anomalistically low BAI) are ecological indicators recording the reaction of the trees to environmental stress factors (Schweingruber et al. 1990). Based on previous studies (Rubio-Cuadrado et al. 2018a, 2018c; Nechita et al. 2019), negative pointer years were determined as those years in which at least 60% of the BAI series of one species showed a decrease in BAI of at least 40% relative to the average BAI in the 4 preceding years. The growth response to these unfavorable years and the ability to recover pre-disturbance growth levels after the disturbance were estimated for the tree-ring series through the resistance (Rt) and resilience (Rs) indices (Lloret et al. 2011):

$$ Rt={BAI}_i/{BAI}_{i- 4} $$

(8)

$$ Rs={BAI}_{i+ 4}/{BAI}_{i- 4} $$

(9)

where BAI_i is the BAI value of the i year. These indices were calculated with 4-year pre- (i-4) and post-disturbance (i + 4) periods using the pointRes package in R (van der Maaten-Theunissen et al. 2015).

2.5 Climate-growth relationships

Since the BAI may contain long-term trends due to non-climate factors, i.e., competition, to determine the main climatic drivers of tree growth, we used ring-width indices (RWI), removing the biological trends of the raw tree-ring widths (TRW) for each species using the ARSTAN program (Holmes 1997). We used a double detrending procedure to remove long-term growth trends (Holmes et al. 1986). First, we fitted a negative exponential curve to the TRW to remove the age effect in the first years of the growth series because of the low competition associated with the open woodland conditions during the establishment period of the analyzed trees, and second, we fitted a 32-year cubic smoothing spline to remove the low-frequency stand dynamics signals. In order to retain climatically related autocorrelation in growth, part of this signal was reintroduced into the residual time series by applying a first-order autoregressive model. Pearson correlation coefficients were then calculated for that period with more than ten samples per species, i.e., for the period 1971–2015, between the RWI (mean series of the detrended TRW) of each tree species and temperature and precipitation grouped on a daily or, alternatively, monthly basis. We tested the Pearson correlations between RWI and all possible climate window widths (periods of year) at daily resolution using the dendroTools package in R (Jevšenak and Levanič 2018). We summarize the climate variables in the selected time window, using mean and sum for temperature and precipitation, respectively. We fixed a 30-day threshold as the minimum widths of the climate window in order to reduce the total number of calculated correlations and therefore the probability of obtaining spurious correlations. Correlation coefficients between RWI and monthly temperature and precipitation were calculated from June of the previous year to November of the growth year considering that the (i) climate during the previous year affects growth during the following year (Fritts 1976) and (ii) growth in the two Quercus species can occur until November (Pérez-de-Lis et al. 2017). Pearson correlations were also calculated between the detrended growth series and the SPEI using time scales from 1 to 20 months (see Section 2.2 “Climate data”).

2.6 Retrospective analyses of species competitive dynamics

To study the specific signals of growth patterns of the three studied species and to analyze their differences, we used the competitive advantage index (CA) based on radial growth data (Rubio-Cuadrado et al. 2018b). This was based on the assumption that for trees of the same species growing in similar conditions of soil, competition, etc., higher growth values imply higher competitive ability (Weber et al. 2008). The CA index was used to compare the percentage of trees of two species showing increasing growth in 5-year intervals. In this way, we are able to study changes in the growth trends (due to climate or other factors) that benefit some species over others. Indeed, a species that maintains a constant downward trend in its relative growth will end up being dominated, benefiting the other competing species. However, this index does not consider other aspects of species competitiveness, such as re-sprouting capacity, carbon storage, seed dispersion capacity, and shade tolerance.

In order to obtain the CA index, firstly, the BAI series of all trees were smoothed, calculating them for 5-year moving intervals since we were interested in the long-term growth trends. Secondly, the number of ascending increments of BAI between consecutive years was summarized for each species. This was then transformed to percentages for species and year (i.e., percentage of trees of F. sylvatica, Q. petraea, and Q. pyrenaica with ascending increments in the 5-year BAI series between two consecutive years) to give a relative series of competitive ability over time. To calculate the competitive advantage series of F. sylvatica over Q. petraea (CA_Fsyl-Qpet), the percentage series of competitive ability for Q. petraea was subtracted from the percentage series of competitive ability for F. sylvatica. In the resulting competitive advantage series, positive values indicate competitive advantage for F. sylvatica, whereas negative values indicate competitive advantage for Q. petraea (Rubio-Cuadrado et al. 2018b). Competitive advantage series for F. sylvatica over Q. pyrenaica (CA_Fsyl-Qpyr) and for Q. petraea over Q. pyrenaica (CA_Qpet-Qpyr) were calculated in the same way.

To calculate the relationships between CA_Fsyl-Qpet, CA_Fsyl-Qpyr, CA_Qpet-Qpyr, and climate, annual series of monthly temperature and precipitation were smoothed using a 5-year moving average. Pearson correlation coefficients between the CA and the smoothed series of climate variables from June of the previous year to November of the growth year were calculated. To calculate the P values of these correlations, we corrected the degrees of freedom, which decrease due to the repeated use of the same data in the moving-average series (Munshi 2016). We also grouped monthly climatic data into seasons to better reflect the relationships between competitive advantage and seasonal climate conditions. The correlations between the CA series and the climatic variables series were analyzed between 1971 and 2015.

2.7 Factors influencing the radial growth of tree species

We used linear mixed-effects models to evaluate BAI trends of F. sylvatica, Q. petraea, and Q. pyrenaica with the following variables included in the saturated models as predictors or fixed factors: DBH of the cored tree at the year of each annual radial growth (annual DBH), as most of the growth variance is explained by tree size (Monserud and Sterba 1996; Cescatti and Piutti 1998; Diaconu et al. 2015); total plot basal area (BA) and the fraction of it belonging to each of the species at the year of each annual radial growth (kriging estimates of annual values at the sampling points), to study the influence of the competition; and the monthly climate variables showing significant correlations with growth for each species (temperatures of March of the current growth year and precipitations from June to October of the previous year and of April of the current year in the F. sylvatica model; temperatures from June to July of the current year and precipitations from July to September of the previous year and from April to July of the current year in the Q. petraea model, and temperatures of November of the previous year and June of the current year and precipitations from July to September of the previous year and from April to July of the current year in the Q. pyrenaica model). Instead of the previous monthly climate variables, we have also tried to introduce, as fixed effects, the temperature and precipitation periods that best correlates with growth for each species using daily resolution (Table 2); however, the resulting models explained less (in the case of F. sylvatica) or similar (in Q. petraea and Q. pyrenaica) variability so to avoid mixing climate data with different temporal resolutions, only the monthly climate data were used in the final models. Only tree identity was considered a random component of the models (1|Tree). We considered all the period with available competition and growth data (1994–2015).

Table 2 Optimal climate windows (window width with the maximum correlation), window lengths in days, Pearson correlation coefficients, and associated probability levels (P) calculated between the ring-width indices (RWI) and climate variables

To predict BAI we adjusted the following linear mixed-effects model with random intercept and fixed slope:

$$ {y}_{ij}=\alpha +{a}_j+{\beta z}_{ij}+{\varepsilon}_{ij} $$

(10)

where y_ij represents log (BAI + 1) for year i and tree j, α is the general intercept, a_j is the random intercept (tree identity), β is the vector of general slopes, z_ij is the vector of fixed effects, and ε_ij is the error with a first-order temporal autocorrelation (AR(1)) structure. We used log-transformed (BAI + 1) data because BAI had a Gamma distribution. The distribution was tested using the Kolmogorov-Smirnov test.

To identify the best-supported model, we constructed all possible combinations of alternative models from the full model considering fixed and random effects and the interactions between the fixed effects. As the restricted maximum likelihood method (REML) estimates an unbiased variance, but does not allow models to be compared by minimizing the Akaike Information Criterion (AIC), we first fitted all possible models using the maximum likelihood method (ML) and then selected the best model by minimizing AIC and finally fitted it again using the restricted maximum likelihood method (REML) (Zuur et al. 2009). The use of the AIC was considered equivalent to using cross-validations to validate mixed-effects models (Fang 2011). The existence of multicollinearity among explanatory variables was evaluated by calculating the variance inflation factor (VIF). VIF values greater than 10 mean that there is high collinearity among variables (Dormann et al. 2013). The percentages of variance explained by fixed (R²_m, marginal R²) and fixed plus random (R²_c, conditional R²) factors were obtained following Nakagawa and Schielzeth (2013).

To study the effect of each variable on tree growth, we used the linear mixed-effects model fitted for each species to predict the evolution of the BAI according to the variation of each variable separately and fixing (as constants) the rest.

3 Results

3.1 Climate trends

The climate data did not show a significant trend when we considered the whole period. However, since 1976, there has been a significant (p < 0.01) temperature increase of + 0.022 °C yr⁻¹ (Fig. 2).

Fig. 2

Climatic trends (annual mean temperature and total precipitation) in the study area (data of the “Pantano El Vado” AEMET station). The slopes of the simple regressions fitted to temperature (from 1976 to 2013) and precipitation (from 1943 to 2014), respectively, are shown

3.2 Growth trends

The three studied species had two clearly differentiated periods according to their growth patterns: an ascending growth trend between 1971 and 1990 and a decreasing (in Q. petraea and more sharply in Q. prenaica) or stable (F. sylvatica) growth trend from 1990 onwards (Fig. 3). The growth changes in F. sylvatica and Q. petraea were very similar over time, their growth curves being almost parallel. However, since 1995, Q. petraea growth has shown a slightly decreasing trend, while F. sylvatica has shown an increasing or stable growth trend, matching the growth of both species between 2007 and 2009, although this suddenly decrease in 2010 and notably since 2015. The spline adjusted to the beech BAI decreases at the end of the period studied due to a sharp drop in growth caused by a spring frost in 2017, which severely damaged the leaves as they were unfolding.

Fig. 3

Growth (BAI, basal area increment) of the three species studied since 1971. We fitted a 30-year-long cubic smoothing spline for each species to highlight BAI decadal to multi-decadal variations. The rigidity of the splines was chosen visually. Note the sharp BAI drop in 2017 observed in beech due to a spring frost

3.3 Climate-growth relationships

The growth of F. sylvatica was directly related to the mean temperature of March of the current growth year (r = 0.29, P < 0.05) and to the total precipitation from June to October of the previous year (r = 0.32, P < 0.05) (Fig. 4). Climate-growth correlations were similar in both Quercus species. Growth of the oak (Q. petraea) was inversely related to the mean temperature from June to July of the current growth year (r = − 0.35, P < 0.05) and directly related to the total precipitation from July to September of the previous year (r = 0.43, P < 0.01) and from April to July of the current year (r = 0.50, P < 0.01). Similarly, the growth of Q. pyrenaica was directly related to the total precipitation from July to September of the previous year (r = 0.33, P < 0.05) and from April to July of the current year (r = 0.44, P < 0.01). Using a daily climate resolution, the periods that showed the greatest correlation with the growth of F. sylvatica were September 29–October 28 for temperature and April 04–May 05 for precipitation (Table 2). In Q. petraea, the periods were June 11–July 17 for temperature and April 04–September 16 for precipitation, and in Q. pyrenaica, the periods were March 17–April 15 for temperature and May 13–September 11 for precipitation.

Fig. 4

Climate-growth relationships based on Pearson coefficients and calculated by relating ring-width indices (RWI) to monthly mean temperature and total precipitation. Months with lowercase letters correspond to the previous year. Horizontal dashed and dotted lines indicate significant correlations at P < 0.05 and P < 0.01, respectively

Only the Quercus species showed significant correlations with the SPEI drought index (Fig. 5). The highest correlations with SPEI corresponded to Q. pyrenaica. The growth of Q. petraea was negatively affected by the cumulative drought severity from March to October of the growth year, whereas the growth of Q. pyrenaica was affected by the cumulative drought severity from April to September.

Fig. 5

Pearson correlations calculated between the ring-width indices (RWI) and the SPEI drought index from January to December (y axes) and for time scales between 1 and 20 months (x axes). Correlation values higher than 0.29 and 0.38 are significant at P < 0.05 and 0.01, respectively

F. sylvatica was the only species which presented pointer years, which were 1995 and 2017. Resistance was 0.45 and 0.18 for the years 1995 and 2017, respectively, and resilience was 0.81 for the year 1995. The resilience for the year 2017 could not be calculated because we only had growth data up to the year 2018. There were spring frosts in both 1995 and 2017.

3.4 Retrospective analyses of species competitive dynamics

During the study period, there were alternating periods of competitive advantage for one or other of the studied species (Fig. 6). In general, we observed no clear trends of the competitive advantage series, except in recent years when the competitive advantage of F. sylvatica over both Quercus species (CA_Fsyl-Qpet and CA_Fsyl-Qpyr) showed a negative trend due to the effects of the spring frost of 2017. However, F. sylvatica showed an increasing growth for a higher percentage of trees than Q. petraea and Q. pyrenaica (i.e., CA_Fsyl-Qpet or CA_Fsyl-Qpyr > 0) in 67% and 65% of the years, respectively, with average values for the competitive advantage series of 5.5% (CA_Fsyl-Qpet) and 10.8% (CA_Fsyl-Qpyr). Although the percentage of trees of Q. petraea with increasing growth was only larger than that of Q. pyrenaica in half of the years, the average value of the CA_Qpet-Qpyr was 5.3%.

Fig. 6

Reconstructed competitive advantage (CA) series of Fagus sylvatica over Quercus petraea (a, CA_Fsyl-Qpet), Fagus sylvatica over Quercus pyrenaica (b, CA_Fsyl-Qpyr), and Quercus petraea over Quercus pyrenaica (c, CA_Qpet-Qpyr) from 1971 to 2018. We fitted a 30-year-long cubic smoothing spline for each CA to show the decadal to multi-decadal variations. The rigidity of the splines was visually chosen

Precipitation from May to June was inversely related with the competitive advantage of F. sylvatica over Q. petraea (Table 3). The temperature from June to July was also inversely related with the competitive advantage of Q. petraea over Q. pyrenaica, while precipitation of February and June was directly related with this advantage. There were no significant correlations of CA_Fsyl-Qpyr with the climate variables.

Table 3 Pearson correlation coefficients and associated probability levels (P) calculated between the different competitive advantage (CA) series and climate variables

3.5 Factors influencing species growth trends

In the three growth models (Table 4 and Fig. 7), the VIF of all predictors was lower than 3, so collinearity was assumed to be low. According to these models, 38.8% (for F. sylvatica), 14.2% (for Q. petraea), and 12.1% (for Q. pyrenaica) of the BAI variability remained unexplained. Similarly, 30.7%, 38.1%, and 36.5% of the variability in each of these species, respectively, were explained by fixed factors (R²_m), and 30.5%, 47.7%, and 51.4% of the variability (R²_c) were explained by differences between trees (soil, orientation, slope, altitude, genetics, vigor, etc.).

Table 4 Regression coefficients of the fixed effects (standard errors are shown in parentheses) and marginal (R²_m) and conditional R² (R²_c) for the best linear mixed-effects models of basal area increment for F. sylvatica, Q. petraea, and Q. pyrenaica, considering the 1994–2015 period

Fig. 7

Effects of interactions between climate variables (P April, April precipitation; T_March, mean March temperature) on radial growth (BAI, basal area increment) of beech according to the selected linear mixed-effects model. Factors are standardized. Shaded bands indicate the 95% confidence intervals. Positions of the data along the abscissa axis are denoted by tick marks

In the beech growth model, there was an interaction between April precipitation and mean March temperature. According to this interaction, when one of the two variables had low values (low temperature or low precipitation), as the value of the other variable increased, growth also increased (Fig. 7). However, when one of the variables showed high values, the other variable had no significant effects (note the width of the confidence interval for this case).

The three BAI models presented a good fit, although in both Quercus models, the fitted BAI showed low inter-annual variability (Fig. 8a). Positive BAI trends of F. sylvatica and Q. petraea were explained by the increase in tree DBH (Fig. 8b). Climate variables explained the inter-annual BAI variability in the three studied species, the BAI decrease in F. sylvatica, and part of the BAI decrease in Q. petraea. Finally, tree-to-tree competition explained most of the BAI decrease in Q. petraea and the BAI decrease observed in Q. pyrenaica.

Fig. 8

Observed and predicted basal area increment (BAI) for Fagus sylvatica, Quercus petraea, and Quercus pyrenaica (a) and BAI prediction series for the variation, within the studied period, of tree diameter (DBH), competition (total BA and the fraction of it belonging to each of the species), and climatic effects separately, keeping the rest of the fixed effects constant (mean value of the period) (b). In the (b) plots, linear regressions were fitted to BAI predictions for the variation in climate variables (keeping DBH and BA constant) to show the long-term BAI trend but not considering annual variability

4 Discussion

Contrary to expectations, we found that the increase in temperatures since the 1970s (Fig. 2) has not benefited the sub-Mediterranean Q. pyrenaica over the temperate species Q. petraea and F. sylvatica. According to our results (Fig. 6), the competitive dynamics benefited Q. petraea and especially F. sylvatica, the growth rates of which are increasing with respect to the other species although Q. petraea currently displays a higher absolute growth (Fig. 3). This result can be partly explained by the increase in competition associated with the cessation of forest management, which has been a major factor behind the decline in growth of Q. petraea and Q. pyrenaica, i.e., they are being outcompeted by beech (Fig. 8b). This study reinforces the notion that the reduction in logging and grazing has been one of the main drivers of forest dynamics in Spain over recent decades (Hernández et al. 2019).

4.1 Climate-growth relationships

Despite the higher vulnerability of beech to extreme climate events such as drought and frost in comparison with both Quercus species (Aranda et al. 1996, 2000, 2005), beech growth was the least sensitive to mean climate conditions (Figs. 4 and 5 and Table 2). Unlike the findings of studies conducted in other temperate forests with Mediterranean influence (Piutti and Cescatti 1997; Rozas 2001) or those of previous studies focusing on aged trees of this forest (Dorado-Liñán et al. 2017a), high temperatures did not seem to be a limiting factor for beech growth. In fact, the only significant negative impact on growth that was identified occurred during the climate window on September 29–October 28 (Fig. 4 and Table 2), when radial growth was negligible or stopped (Martinez del Castillo et al. 2018). Similarly, site conditions did not seem to be especially limiting to beech growth. While many missing rings were detected in dominant beeches in other temperate forests located further north due to adverse climate conditions or suppression caused by intense competition (Rubio-Cuadrado et al. 2018c), no missing rings were detected in this case. These results agree with those of other recent studies which point to adaptation to drought in the most southerly beech forests (Tegel et al. 2014; Cavin and Jump 2017; Muffler et al. 2020).

The influence of precipitation from June to October of the previous year on beech growth (see Section 3.3 “Climate-growth relationships”) is related to the amount of carbohydrates stored in the previous growing season and therefore available at the beginning of the current vegetative period (Bréda et al. 2006), whereas the influence of March temperature could be related to earlier breaking of winter dormancy along with enhanced stimulation of the cambium activity (Seidling et al. 2012). In this respect, beeches in the study area sprout at the end of April (Millerón et al. 2012). It is also worthy of note that the beech trees in this study were not sensitive to June precipitation, which has been shown to drive growth in this species in populations of northern Spain (Rozas et al. 2015), while the precipitation just before and at the beginning of the vegetative period was significant (Table 2). This may be explained by the fact that the phenology of these trees was more advanced in comparison with northern populations located in cool, wet sites. Finally, spring frosts cause noticeable growth losses in such populations (Dittmar et al. 2006). In the study area, four spring frosts have occurred in the last 25 years: 1995, 2010, 2013, and 2017 (data not shown). The years 1995 and 2017 were the pointer years, but sharp drops in the beech radial growth was found in all these years.

Radial growth of the two Quercus species was also related with the precipitation during the previous growing season, as well as with spring and summer precipitation of the growth year. In addition, the temperatures of June and July (see Section 3.3 “Climate-growth relationships”), when the highest radial growth occurs (assessed through band dendrometer measurements not shown in this paper), and those from mid-March to mid-April (Table 2), just before the beginning of the vegetative period, were inversely related with the total growth of Q. petraea and Q. pyrenaica, respectively, which may be due to their higher proportion of living parenchyma in sapwood and higher respiration values per unit volume of living parenchyma compared with F. sylvatica (Rodríguez-Calcerrada et al. 2015).

4.2 Species competitive dynamics

The BAI trends (Fig. 3) and the competitive advantage series (Fig. 6) are clearly linked: the higher growth increments (i.e., the higher slope of the splines in Fig. 3) of F. sylvatica and Q. petraea compared with Q. pyrenaica until 1990 led to large positive values, around 20%, of CA_Fsyl-Qpyr and CA_Qpet-Qpyr. From then on, there was a stabilization of F. sylvatica growth along with a decrease in growth in both Quercus species, and consequently, CA_Fsyl-Qpet and CA_Fsyl-Qpyr continued to maintain positive values, while CA_Qpet-Qpyr alternated between positive and negative values, although the positive values predominated due to the greater decreases in Q. pyrenaica growth. Finally, there was a sharp drop in CA_Fsyl-Qpet and CA_Fsyl-Qpyr due to damage suffered by beech in the late frost of 2017.

Most of the relationships between the competitive advantage series and climate (Table 3) can be explained by the different climate-growth relationships presented in Fig. 4, which are associated with differences in spatio-temporal use of water, cambium phenology (Michelot et al. 2012), and physiology (Aranda et al. 1996, 2005). While beech growth showed no relationship with precipitation from May to June, both Quercus species presented significant positive relationships. Thus, the precipitation in these months benefits the Quercus species over F. sylvatica. Similarly, Q. pyrenaica growth showed no relationship with summer temperatures, whereas Q. petraea presented significant inverse relationships. Consequently, high summer temperatures, along with low precipitation in February and June, will have a negative impact on Q. petraea against Q. pyrenaica. High precipitation in February may benefit Q. petraea growth to a greater extent than Q. pyrenaica because its radial growth starts earlier, as in similar Iberian temperate forests (Pérez-de-Lis et al. 2017),

In our hypothesis, we expected an improvement in the competitive capacity of the more drought resistant species, as previously reported in temperate forests of northern Spain under mesic conditions (Rubio-Cuadrado et al. 2018b). However, in recent decades, the competitive advantage series has shown a tendency favoring temperate species, especially beech, despite the negative impact on beech growth of extreme climate events such as droughts (Rubio-Cuadrado et al. 2018c) and spring frosts (Dittmar et al. 2006). These results may be related to the lower sensitivity and greater resistance to drought shown by rear-edge beech populations compared with those situated at the core of the species range (Cavin and Jump 2017; Serra-Maluquer et al. 2019). This greater resistance may be explained by favorable local environmental factors (e.g., site topography, soil water holding capacity) or by acclimation to drought (Rodriguez-Calcerrada et al. 2010) through phenotypic plasticity or genotypic adaptation. In any case, recent studies have revealed high resilience to drought in beech saplings, which rapidly resume physiological activity after the dry period (Pflug et al. 2018). These aspects clearly merit further investigation to determine whether adult beech trees display similar resilience to different extreme climate events (droughts, frosts) in rear-edge populations.

4.3 Factors influencing species growth trends

The model has been built for the period 1994–2015 as we do not know the precise evolution of the competition before the first forest inventory carried out in 1994. Although the time frame considered may be short, due to the relatively intensive use by cattle in the past and its subsequent abandonment, there has been a considerable increase in competition in this short period of time, which makes this forest an ideal site to study its effect on growth. The systematic sampling of trees performed, where all plots (of the inventory carried out in the studied forest in 2005) with mixed forest were sampled, and the high EPS values (Table 1) make our models representative of the growth dynamics of the dominant and co-dominant trees of the whole forest. However the short time frame considered together with low number of trees per species sampled limits the scope of the results. Further research studying longer periods of time and larger number of trees, and covering the distribution area of the species, are needed to know the factors that are influencing at a global level the growth trends of F. sylvatica, Q. petraea, and Q. pyrenaica.

As reported in other studies (Monserud and Sterba 1996; Cescatti and Piutti 1998; Diaconu et al. 2015), most or at least a large proportion of the tree growth in F. sylvatica and Q. petraea was explained by tree size (DBH; Table 4). In alternative models for Q. pyrenaica that included the DBH (results not shown), this variable also explained most of the radial growth variability, but was not significant because it had a very wide confidence interval. Both Quercus BAI models, especially the Q. pyrenaica model, were less sensitive to climate than the F. sylvatica model, as shown by the lower coefficients of the climate factors (Table 4), hence displayed lower inter-annual variability (Fig. 8a).

According to our BAI models (Table 4), the increase in competition following the cessation of forest management is the most important factor related with the negative growth trend of Q. petraea and Q. pyrenaica, although climate also explained some of the Q. petraea growth decline (Figs. 3 and 8b). These results are concordant with previous works, where high-frequency growth changes are related to climate, while radial growth is explained by the size and competition (Sánchez-Salguero et al. 2015; Liang et al. 2019). Indeed, competition among trees may cause much larger reductions in forest growth and carbon uptake than climate stress (Vayreda et al. 2012; Rozas 2014). The greater importance of competition rather than climate in explaining tree radial growth has also been reported in large-scale studies conducted both in temperate-boreal (Foster et al. 2016) and in water-limited ecosystems (Gómez-Aparicio et al. 2011). However, some studies give greater weight to the climate and global warming effect on growth trends, especially in Mediterranean rear-edge forests (Piovesan et al. 2008; Dorado-Liñán et al. 2017b; Charru et al. 2017; Dorado-Liñán et al. 2019), although these studies do not usually analyze the effect of the evolution of competition because it is uncommon that inventory data is available for a sufficiently long period of time. Our study shows the greater relative importance of competition in explaining radial growth not only in Mediterranean Q. pyrenaica but also in a rear-edge population of Q. petraea, providing further evidence of the complex relationships between climate and growth decline in some temperate forests at the equatorial range edge of their distribution (Tegel et al. 2014; Cavin and Jump 2017; Muffler et al. 2020). In fact, beech is currently expanding in Spain despite the reported aridification trends (Sánchez de Dios et al. 2016, 2020) due to its competitiveness, shade tolerance, and greater dispersion capacity (Harmer 1994; Millerón et al. 2013), which allows it to access the best situations in terms of soil depth and orientation. The competition was not significant in the F. sylvatica BAI model (Table 4), possibly due to its greater competitiveness in comparison with the oak species (Hein and Dhôte 2006; Manso et al. 2015). This could result in a negative influence of beech admixture on oak. However, according to the models obtained, the presence of Q. pyrenaica seems to be the most limiting for the growth of both Quercus, while the presence of Q. petraea was the least limiting. This result may be influenced by the larger average size of Q. pyrenaica trees (see Section 2.1 “Study area”) and the direct relationship between size and competition (Stadt et al. 2007). Furthermore, previous studies have shown that competition for resources in unfavorable climatic years may be lessened in mixed beech-oak stands with respect to pure stands (del Río et al. 2014).

Changes in tree-to-tree competition must be considered in order to understand long-term changes in growth and forest dynamics since they may be more relevant than climate to forecast changes in forest composition and dynamics (Sánchez-Salguero et al. 2015). In addition, changes in competition can drive or modulate the responses of tree growth and functioning (e.g., water-use efficiency) to climate (Linares et al. 2010; Martín-Benito et al. 2011; Fernández-de-Uña et al. 2015; González de Andrés et al. 2018). The effect on tree growth of the increase in competition resulting from the abandonment or change in previous forest management practices (which is common in many European forests) and that of global warming can be difficult to unravel when analyzing the causes of growth decline. Therefore, by combining dendrochronological data with historical forest inventories, it is possible to provide estimates on the evolution of stand competition among trees over time and improve our knowledge on the causes of growth decline in some of the temperate rear-edge forests.

5 Conclusions and management implications

There has been an increase in temperatures in the study area from 1975 onwards. However, the rear-edge Q. petraea and especially F. sylvatica presented low sensitivity to climate and showed favorable evolution of their competitive capacity over the Mediterranean Q. pyrenaica. F. sylvatica also showed higher competitive capacity than Q. petraea throughout most of the studied period. Spring frosts, in the case of F. sylvatica, and increased competition resulting from the cessation of forest management, for the two Quercus species, were the most limiting factors to tree growth. Management strategies, such as selective thinning and promoting of sexual reproduction in Q. pyrenaica, should be considered a means to improve long-term growth and enhance the resilience of Quercus populations.