The study area is located in Aspurz (42º 42’ 31’’ N, 1º 08’ 40’’ W), in the western Spanish Pyrenees, Navarre (Figure 1). It is located in one of the most productive Scots pine forests in Spain (Puertas, 2001), which site index is 29 m at 80 years. The forest is an even-aged mixed stand dominated by Scots pine, naturally regenerated in the mid-1960s after strip-like clear- cutting, with mean dominant tree height of 20.4 m in 2014. Beech is the second most abundant tree species, being very abundant in the understorey (< 2 m) at certain mixed areas. However, it also occurs as an upper continuous canopy layer including the four crown classes (Smith, 1986), where competition between beech and pine crowns is evident. Additionally, some isolated beech individuals can be found between this upper layer and the understorey. In 2008, beech had a mean crown cover of ca. 38% and mean beech age was 35-40 years. The study site is North-oriented and it is located at a mean altitude of 642 m, a mean slope of 7 % and it develops on a Haplic Alisol (Blanco, 2004). The climate is classified as a cold wet Mediterranean type with water deficit usually in July and August, and frequent frosts in winter and early spring. Mean annual precipitation is 921 mm and mean annual temperature, 12.0 ºC (data from nearby Navascués station located at 2.7 km from the study site, 42º 43’ 06’’ N, 1º 06’ 55’’ W, 615 m; period 1984- 2014).

Nine rectangular plots (30 m x 40 m) were included in the experimental design, installed by the Servicio Forestal del Gobierno de Navarra following a complete randomized block design (Andrew, 1986). Three different thinning intensities (3 thinning intensities x 3 replicates) were applied in November 1999 (0, 20 and 30% basal area removed, hereafter unthinned-U -, moderate-M- and severe-S- thinning) and March 2009 (0, 20 and 40%). The first thinning was carried out removing mainly suppressed or intermediate trees, and some dominant or codominant trees with malformed stems. During the second thinning, mainly subdominant or dominant trees were removed, as most trees were in the upper strata as the result of the stand dynamics and the first thinning. The highest thinning intensity was increased up to the 40 % of the basal area during the second thinning in accordance to the silvicultural trends applied on similar surrounding stands at that moment. Similarly, following the silvicultural guidelines applied in Navarre on mixed stands with beech as secondary species, only Scots pine trees were thinned. In 1999, stumps were left on the ground, and logs and most branches from the felled trees were left outside the plot limits. In 2009, stumps and branches, the latter cut into pieces, were left on the forest floor but the stems, including the bark, were removed from the experimental area. The same treatments were applied in a buffer zone of 5-10 m around the plots.

Within each plot, two discontinuous subplots can be distinguished based on the crown cover of beeches taller than 2 m: mixed beech- pine and pure pine subplots (hereafter abbreviated as MC and PC subplots). Six treatments are thus determined: UM (unmanaged mixed beech-pine), UP (unmanaged pure pine), MM (mixed beech-pine in moderate thinning), MP (pure pine in moderate thinning), SM (mixed beech-pine in severe thinning) and SP (pure pine in severe thinning). The experimental design can be thus specified as a split-plot (von Ende, 2001) with three replicates of six subplots (Table 1). Data on beech tree cover is not available for 1999, but comparison between beech cover data in 2003 and 2008 suggests that mixed canopy area has remarkably spread over the plots since the onset of the experiment.

The study site was considered as a pure pine forest due to the low presence of beech by the foresters who set up the experimental plots in 1999. However during the course of our studies in the site we have shown that beech cover is expanding in the stand and that beech has a great effect on forest nutrient cycling (see Primicia et al., 2014). Therefore, the main aim of this study was to study the effect of beech presence on the dynamics of Scots pine growth. Pure beech plots have been lately installed near the experimental site to model the development of whole stand growth in future analyses.

In 1999, all alive and standing dead Scots pine trees within the plots were labelled and their relative coordinates calculated using a compass and the Vertex III hypsometer (Haglöf, Sweden). Two perpendicular measurements of their diameter at breast height were taken and the measuring height was marked to reduce the error when comparing with the following measurements. The inventory was repeated in 2009 (before the second thinning) and in 2014. The variables measured during each inventory are shown in Table 2. Stem volume equations were derived from a sample of 229 Scots pine trees covering all diameter classes felled down in 1999 and 2009 (n = 200 and n = 29 trees for 1999 and 2009, respectively). The best-fitted model was selected as that showing the lowest Akaike Information Criterion (Burnham & Anderson, 2002) and highest R ² adj values once the assumptions of normality and homoscedasticity were satisfied (Suppl. Table S1 [online supplement]). To analyse Scots pine mortality in the study plots, we used data on standing dead trees (snags) recorded in 1999 and data on both standing and lying dead (logs) trees recorded in the 2009 and 2014 inventories. Mean mortality rates (frequency, basal area) and periodic annual increments (PAI) in basal area (PAIBA), stand volume (PAIV), and in tree volume (iv̄) were calculated for the two periods determined by the three inventories (1999-2009 and 2009-2014).

All statistical analyses were performed using the R software (R Core Team, 2013). The effects of thinning intensity and type of canopy on the tree and stand structure variables (density, N; basal area, BA; stand volume, V; tree volume, v̄ data before thinnings), PAIBA, PAIV, iv̄, and on mean mortality rates, were assessed using linear mixed effects models following a split-plot design with a repeated measure (Pinheiro & Bates 2000) using the nlme package (Pinheiro et al., 2009). We used thinning intensity, canopy type and the inventory year (or period for PAI variables and mortality rates) as fixed effects. To analyse the PAI variables, we also included the quadratic mean diameter (Dg) into the models as a fixed effect, to account for the influence of the stand structure on growth and mortality rates. We choose Dg because it was the stand variable most correlated with the PAI variables. We did not include the dominant height, like in other studies, because due to the proximity and situation of the study plots on the same slope, no differences in site quality can be expected. As random effects, we included subplot (for variables at the stand scale) or tree nested in subplot (tree variables), since they can be considered as a random representation of the population. An autoregressive correlation structure of first order was included to account for the repeated measures on the same subplot or tree (Pinheiro & Bates, 2000). For instance, to analyse iv̄, a variable calculated at the tree scale, the proposed model was stated as follows:

iv̄ = I + Dg + thinning * canopy * period + tree:subplot + corAR (tree:subplot)

where: iv̄ represents the Tree volume annual increment; I is the intercept; Dg is the subplot quadratic mean diameter; thinning and canopy represent the study treatments; period is the sampled period for which iv̄ has been estimated; tree:subplot, a random term specifying the effect of individual trees nested within each subplot; and corAR (tree:subplot), the autoregressive correlation structure to account for the repeated measures on each tree within each subplot.

The appropriateness of the random and correlation structures was analysed by comparing nested models with and without the structures with the likelihood ratio test using the restricted maximum likelihood estimation procedure. The significance of the fixed effects was analysed by comparing nested models with and without the fixed terms using the maximum likelihood estimation procedure (Zuur et al., 2009).

To analyse whether the effects of treatments on the Tree volume annual increment (iv̄) depended on tree size, we used linear mixed effects models fitted to the data set split into four classes based on the 25%, 50%, 75% and 100% percent quartiles (Q₂₅, Q₅₀, Q₇₅, and Q₁₀₀, respectively) of the tree diameter at breast height (DBH) at the beginning of the study. Thus, the data set was split into four classes: Q₂₅, Q₅₀ - Q₂₅, Q₇₅ - Q₅₀, and Q₁₀₀ - Q₇₅. The proposed models were similar to the one stated above, but in this case we used the initial DBH instead of quadratic mean diameter as covariate.

Evolution of the main stand variables in the different thinning treatments are shown in Suppl. Table S2 ([online supplement]). Scots pine density and basal area were lower in thinned than in unthinned plots only after the second thinning (P < 0.001, Table 3, Figure 2a), while stand volume was higher (P <0.01, Table 3) in U than in M and S plots, although differences were mainly observed in 2014 (295 ± 14, 253 ± 17 and 222 ± 11 m³ ha^-1 in U, M and S plots, respectively, Figure 2a). Tree volume decreased with increasing thinning intensity at the beginning of the study, but in 2014, it was highest in the M plots (P < 0.001, Table 3, Figure 2a). Scots pine basal area and stand volume was significantly (P <0.001) higher in pine than in mixed canopy (Figure 2b). Differences between canopy types for pine density decreased over time (PC > MC, P = 0.09, Figure 2b), except for tree volume, which was similar between canopy types at the beginning of the experiment, but it was higher in PC than in MC at the end of the study (P < 0.05, Figure 2b).

In 1999, the number of snags was similar in the severe than in the unthinned plots (Figure 3). Mortality rates (frequency) was significantly influenced by the interaction between thinning intensity and period (P < 0.05, Table 3), since higher rates occurred in U plots compared to M and S plots during the first period, but higher values were observed in U and M plots than in S plots during the second period (Figure 3). Basal area mortality rates increased from S plots, followed by M plots, and reaching the highest values in U plots (P < 0.001).

Both periodic annual increments in basal area and in stand volume decreased with increasing the quadratic mean diameter (P < 0.001, Table 3), and they were significantly higher in the thinned than in unthinned plots (P < 0.001). In unthinned plots, a decrease in the stand volume occurred after the second thinning due to natural pine mortality.

Tree volume annual increment (iv̄) was always lowest in U plots, and while it was similar in thinned plots during the first study period, it was lower in S than in M plots during the second study period (P < 0.001, Table 3, Figure 4a). Tree volume annual increment was similar between canopy types during the first study period, but it was lower in MC than in PC during the second period (Figure 4b). Nevertheless, the reduction of iv̄ in the MC was not observed in the thinned plots (P =0.07, Table 3, Figure 4c). The effect of treatments on iv̄ was influenced by tree size, being the biggest trees (DBH comprised between Q₅₀ and Q₁₀₀) the most influenced ones by treatments (Table 4, Figure 5). Tree volume annual increment of the smallest trees (Q₂₅) showed practically no differences among thinning intensities, while bigger trees (DBH comprised between Q₅₀ and Q₁₀₀) showed higher iv̄ in M plots, followed by S plots, and finally, by U plots, especially during the second study period (Figure 5). Tree volume annual increment of all trees except the smallest ones was higher in PC than in MC, although the effects frequently depended on the study period, as differences between canopy types increased after the second thinning.

The tree volume annual increment (iv̄) of Scots pine was slightly enhanced after the first thinning, but no differences were observed between M and S plots. Primicia et al. (2013) also found a weak and temporary effect of the first thinning on the tree ring width of dominant and codominant pines in the same study site. Although thinning from below in Scots pine stands appears to particularly enhance the growth of intermediate trees (Río et al., 2008), we observed that the biggest trees (DBH comprised between Q₅₀ and Q₁₀₀) were the most responsive ones to the first thinning. In the study site, the initial stand structure of Scots pine was greatly variable among plots (lower pine density with higher quadratic mean diameter, higher tree and stand volume in unthinned than in thinned plots). Additionally, high natural pine mortality occurred in U and M plots after the first thinning, highlighting a strong competition among trees. Both the initial variability in the stand structure and the disparities in the mortality rates may have disguised to some extent the effects of the first thinning on iv̄.

After the second thinning, iv̄ was highest in the M plots, especially in the case of the biggest trees, contrastingly to the results obtained in previous studies, where increment in diameter increased with increasing thinning intensity (Montero et al., 2001; Mäkinen & Isomäki, 2004). In Mediterranean sites, thinning could enhance tree growth by improving the tree water status (Bréda et al., 1995; Ma et al., 2010; Molina & del Campo, 2012), since tree growth is usually mainly limited by water deficit (Sabaté et al., 2002). However, in our study plots, soil moisture was higher in unthinned than in thinned plots even 9 years after the first thinning, probably due to higher transpiration rates in the latter (Primicia et al., 2013). Additionally, wind speed within the stand and evaporation may also increase following thinning, enhancing the transpiration rates of individual trees (Aussenac, 2000). Therefore, the trees water balance might have been negatively affected by the severe- thinning treatment. The higher iv̄ in M than in S plots after the second thinning could be related to a fertility decrease because of nutrient loss in the highest thinning intensity. In the study site, the main nutrient limiting Scots pine growth appears to be phosphorus (Primicia et al., 2014), whose reserves are susceptible to overexploitation (Blanco et al., 2005). In this context, we suggest that the decrease in soil fertility is mainly due to a reduction in litterfall after thinning. For instance, during the first 30 months after thinning in 1999, the percentage of P returns for total aboveground litterfall in thinned relative to control plots ranged from 16.9% (M plots) to 22.3% (S plots; Blanco et al., 2008). Prescriptions for forest thinning in high-productive Scots pine stands following the economic criteria and adapted to our region (exploitation cycle of 80 years, first thinning at 20- 25 years, rotation length of 10 years, removal of 30% of the basal area; Río & Montero, 2001; Puertas, 2003) could indeed result in a significant net phosphorus loss in the study site at the long term (Blanco et al., 2005). Current forest thinning practices applied on surrounding similar stands should be consequently verified from an ecological perspective at the long term.

Thinning treatments enhanced the periodic stand volume annual increments, even though no differences between M and S plots were found. The differences between unthinned and thinned plots were especially notable after the second thinning, as a decrease in the stand volume occurred because of natural pine mortality in the unthinned plots. In the thinned plots, the yield reduction by thinning is counteracted by the enhancement of individual tree growth, resulting in lower stand density but bigger tree size. However, a reduction in either stand or merchantable volume and stand growth has been frequently observed following moderate to very heavy thinning, although it was also absent after light thinning or in control plots in other Scots pine stands few years after thinning (Montero et al., 2001; Mäkinen & Isomäki, 2004; Mäkinen et al., 2005; Río et al., 2008; Crecente-Campo et al., 2009).

Beech presence apparently caused a reduction in the Tree volume annual increment, agreeing with previous results (Primicia et al., 2013; Río et al., 2014). In our study site, characterized by summer water deficit, the negative effect of beech on iv̄ may have been related to the beech root system, since it is highly competitive for water and nutrient resources (Curt & Prevosto, 2003). The importance of light in competition in mixed Scots pine-beech forests in the same region has been also highlighted (Río et al., 2014), where size-asymmetric competition (resource uptake depends on the relative sizes of the competitors) is apparently stronger than size-symmetric competition (resource uptake is independent of the relative sizes of the competitors, Schwinning & Weiner, 1998). The reduction of the Tree volume annual increment due to beech presence gradually increased during the study period, related to beech crown development. Nevertheless, this decrease in iv̄ in the mixed compared to the pure canopy did not happen in the thinned plots. Primicia et al. (2013) also observed that Scots pine growth was higher in PC than in MC in the same study site 9 years after the first thinning, but only in the unthinned plots. They suggested that these results were related to differences among treatments (thinning intensity x canopy type) in the tree water status, tree size and tree-to-tree competition.

The differences in Scots pine density between canopy types found at the onset of the experiment (PC > MC) decreased over time. This pattern could have been partly due to Scots pine mortality in the MC, or due to an unequal thinning intensity within plots. Although the thinned trees were broadly randomly distributed within the plots, some Scots pines near the beech trees were maintained to better study the inter-specific processes between both species and to homogenize the stand structure between canopy types.