Mixing effects on growth efficiency in mixed pine forests

ForestryAn International Journal of Forest Research Forestry 2017; 90, 381–392, doi:10.1093/forestry/cpw056 Advance Access publication 9 November 2016 Mixing effects on growth efficiency in mixed pine forests José Riofrío 1 *, Miren del Río 1,2 and Felipe Bravo 1,3 1 Sustainable Forest Management Research Institute, University of Valladolid & INIA, Av. Madrid 44, Palencia 34004, Spain 2 Department of Silviculture and Forest Management, INIA, Forest Research Centre, Ctra. A Coruña, km 7.5, Madrid 28040, Spain 3 Departamento de Producción Vegetal y Recursos Forestales, E.T.S de Ingenierías Agrarias, Universidad de Valladolid, Palencia, Spain *Corresponding author. E-mail: [email protected]s Received 4 April 2016 Increased interest in mixed forests is due to evidence of them being more resource-use efficient and stable forest systems. However, intrinsic and extrinsic factors moderate interspecific species interactions generating different effects in productivity. Here, we explore a method to detect mixing effects in a specific mixture combination (Pinus sylvestris L. and Pinus pinaster Ait.), comparing the growth of mixed stands with that of monocultures. Combined tree and stand-level analyses also helped determine which mixing effects are most important for forest functioning and how changes at one level influence patterns at another level. Data from the Spanish National Forest Inventory were used to compare growth efficiency in mixed and pure stands; we relied on relative stand density indices to determine species-specific site occupancy. This same concept was used to evaluate competition status and inter/intra-specific competition effects as modifiers of potential growth at the tree-level. We observed that growth efficiency in both species increased with the proportion of the complementary species in the stand. At the tree-level, intraspecific competition was higher than interspecific competition in Scots pine tree growth, showing that it had benefited from the mixture. In contrast, maritime pine did not show a competitive response to the interspecific interaction, indicating that tree growth was more strongly influenced by the competition structure (size-symmetric and size-asymmetric) than by the species of the competitors. Our results highlight the importance of combining stand-level analysis with that of tree-specific competition relationships when studying mixed-species forests. Introduction The fact that mixed forests can be more productive than single species stands has been observed in many species combinations and ecosystems (Kelty, 2006;Vilà et al., 2013;Forrester, 2014). However, the generalization that forest productivity increases with species mixture is inaccurate, since interspecific species interactions can be moderated by site conditions (Huber et al., 2014), climate (Manso et al., 2015), stand density (Condés et al., 2013) or species assemblage (Toïgo et al., 2015), which affect productivity (Pretzsch and Schütze, 2016), mortality (Condés and Río, 2015) or resource-use efficiency (Richards et al., 2010; Forrester, 2014). Interactions between any two species can also change along spatial and temporal gradients, adding further complexity (Forrester, 2014). Most studies of mixed forest in temperate zones focus on mixtures that combine species with different traits (shade tolerance, height growth rate, crown structure, leaf phenology, root depth) (Kelty, 2006;Forrester, 2014), due to the higher probability of complementarity (reduced competition and facilitation) in the allocation and efficient use of available resources between species. Because competition between species increases as species niches overlap, when organisms occupy the same spatial and temporal habitats and use the same resources, the strength of interspecific competition is affected and they could interact as the same species (functional redundancy concept) (Naeem, 1998;Rosenfeld, 2002). However, a positive diversity-productivity relationship has been also observed in discrete species functional groups or homogeneous forest types, such as alpine coniferous or Mediterranean conifer forest, showing higher wood production in the mixed-species stands than in the monospecific stands (Vilà et al.,2013 ;Jucker et al.,2014). Nevertheless the ecological mechanisms underlying the positive relationship between tree species richness and wood production remain unclear. Along these lines, combining tree and stand-level analyses helps to determine which mixing effects are most important for forest functioning and to examine how changes at one level influence patterns at another level (Forrester and Pretzsch, 2015). While tree-level analysis can be directly carried out through regression techniques where the response is single-tree growth, mixing effect at the stand-level is usually analyzed through the possible growth deviations in mixtures as compare to neighbouring pure stands. In this respect, species proportion was used as a proxy for the partitioning of growing space: the allocation of resources such as light, water and nutrients among the species in the stand (Río et al., 2016). ©Institute of Chartered Foresters, 2016. All rights reserved. For Permissions, please e-mail: [email protected]. 381 Downloaded from https://academic.oup.com/forestry/article/90/3/381/2452295 by UNIVERSIDAD DE VALLADOLID user on 23 November 2021 The species proportion may vary considerably depending on the species proportion definition, leading to varying net effects on stand growth (Pretzsch, 2009;Huber et al., 2014). Selecting a species proportion definition that allows for different potential density between species might therefore be more advisable, since it provides a better estimate of the growing space occupied by each species (Sterba et al., 2014). This approach has been used to analyse mixture effects at stand and tree level based on forest inventory data (Río and Sterba, 2009;Condés et al., 2013;Huber et al., 2014;Río et al., 2014;Condés and Río, 2015) because some factors influencing growth are not well described and direct comparison among pure and mixed species plots is not possible. When large differences are observed in species potential density, relative species proportions might provide more reliable mixing effects and more accurately reflect the net total mixing effect as well as the relative importance of under/ overyielding by species (Sterba et al., 2014). However, this approach requires potential density estimates for each species. Competition from neighbours is an important driver of tree growth and might cause much larger reductions in potential growth than climate or tree size (Sánchez-Salguero et al., 2015). Moreover, sensitivity to competition (competitive response) in target trees varied markedly along environmental gradients, causing significant rank reversals in species performance, particularly under xeric conditions (Gómez-Aparicio et al., 2011). Scots pine (Pinus sylvestris L.) and maritime pine (Pinus pinaster Ait.) are two of the main forest species in Spain (1.20 and 0.68 million ha, respectively) and grow in pure and mixed stands, either naturally or derived from species selection for afforestation (Serrada et al., 2008). In addition to their wide distribution and forest area, they hold great ecological and socio-economic value. Mixed stands where these two species coexist are particularly interesting because of their location at the edges of range distribution of P. sylvestris, where ecological conditions (high temperatures, frequent droughts) approach the species tolerance limit and the most drastic effects of climate change are predicted (Matías and Jump, 2012). The aim of this study was to combine tree- and stand-level analysis to evaluate the mixture effect in Scots pine and maritime pine mixed stands, in an attempt to understand the species-specific interactions that favour higher profitability in the mixture than in monospecific stands. For this purpose we established the following approach: (i) Determine the species-specific potential density and space occupancy. (ii) Evaluate if the growth efficiency at the stand level is modified by mixture in terms of the growing space occupied. (iii) Test different competition structures to evaluate the influence of intraspecific and interspecific competition on tree growth. Materials and methods Data Data from the Second and Third Spanish National Inventories (NFI) were used to determine the areas where pure and mixed stands of maritime pine (P. pinaster) and Scots pine (P. sylvestris) are aggregated (Figure 1). NFI plots consist of four circular concentric subplots with radii of 5, 10, 15 and 25 meters, where diameters and heights of all trees over 7.5, 12.5, 22.5 and 42.5 cm breast height diameter, respectively, are recorded. Estimated stand variables (per hectare) included the number of stems (N), basal area (G), dominant height (H O ), quadratic mean diameter (dg) and total stand volume (V) in pure and mixed plots; were calculated using adequate expansion factors. Individual-tree stem volume over bark was estimated from volume equations proposed by the Second NFI (ICONA, 1990). The origin of the stands, natural (N) or plantation (P) was also included; maps with stand origin data are available at https://sites.google.com/site/sigforestspecies/home. Plots from both inventories were selected based on species composition. Pure plots of Scots pine and maritime pine stands were considered only when Gof the target species was higher than 90% of the total. From these plots, mortality rates between inventories that indicated self-thinning conditions, were used to obtained the dataset to fit the self-thinning boundary line (Reineke, 1933) (Supplementary data, Table S1). Thus, we estimated the maximum stand density index value (SDI Max ) for each species, as described in the ‘Stocking proportion and maximum stand density’section. Both species were present in the mixed plots selected, and the combined proportion of Gfor both species accounted for at least 90% of the total. The proportion of other species remained lower than 10%, and the proportion of each target species in the mixture plots was higher than 15%. With these criteria, we sought to avoid effects from other species in the mixed plots and consider the full range of mixture proportion (Supplementary data, Figure S1). Only plots measured in both inventories were considered to analyze the effect of the mixture on growth efficiency. Plots with more than 5% of the total basal area and by species removed between the two consecutive inventories were excluded. Volumen increment was calculated by the difference in standing volumen between inventories, extracting natural mortality and including in-growth from the Third NFI. Plots lacking re-measurement, origin data, with damages due to biotic or abiotic factors (G>40% affected) or with recorded growth anomalies were also rejected. Finally, to reduce the effects of different abiotic growing conditions between mixed and pure plots: only pure plots within a buffer of <5km from the mixed plots were included. To avoid the use of very low density plots, those with SDI below the 15% SDI max were excluded in the subsequent analysis. A total of 115 mixed plots, 210 pure Scots pine plots and 148 pure maritime pine plots were selected and used to analyze growth efficiency (Table 1). Climatic data (mean monthly, seasonal and annual temperature and precipitation values from 1951 to 1999) were added using available raster maps (at 1 ×1 km scale) from a functional phytoclimatic model (Gonzalo Jiménez, 2010). The corresponding phytoclimatic regions were assigned to each selected plot according to their coordinates. These variables were used to test if the main environmental conditions influenced the mixture effects in the stand model. Stocking proportion and maximum stand density To compare growth in mixed and pure stands, negative or positive interaction effects on growth must be related to the area occupied by the respective species. The area available for a species in a mixed stand (species proportion by area) can be defined as the ratio of the observed space occupancy and the potential (maximum) space occupancy for this species and site. The maximum stand density index (Reineke, 1933) and the maximum basal area (Sterba, 1987) have well-developed theories that describe potential density and can be estimated from large inventory data (Río and Sterba, 2009;Condés et al., 2013). They are expressions of the same ecological principal (Vospernik and Sterba, 2015) and both approaches produced similar estimations of mixture effects when the potential densities of the species did not differ by Forestry 382 Downloaded from https://academic.oup.com/forestry/article/90/3/381/2452295 by UNIVERSIDAD DE VALLADOLID user on 23 November 2021 much (Sterba et al., 2014). In this way, the ‘stocking proportion’or species proportion by area (P i ) (equation 1) can be estimated using the respective relative stand densities of the target species (SDIR i ) (equation 2)(Sterba et al., 2014): =+()PSDIR SDIR SDIR 1 ii ij () == ⋅() N C SDIR SDI SDI 25 2 iiidg E iE max 25 i i i where P i is the proportion of species iby area, SDIR i is the observed relative stand density of species i , SDI i is the stand density index of species i and index jis referred to competitor species. The denominator in equation (2) represents the maximum or potential stand density index ( S DIma x ) of species i in a fully stocked pure stand. Thus, mixing proportion were adjusted to take account interspecific differences in carrying capacity at a given site and allowing a plausible comparison of productivity per hectare between pure and mixed stands (Forrester and Pretzsch, 2015;Río et al., 2016). In order to use this approach to calculate the stocking proportion of a given species, the maximum stand density had to be determined for both species. To derive the potential density of the species in pure stands we used the maximum stand density index S DIma x as the self-thinning rule (Reineke, 1933), so C and E in equation (3) depend on the species: =⋅ ()NCdg 3 E max The parameters of the self-thinning boundary (equation 3) for both species were estimated with a log-linear quantile regression (QR) in equation (4) using the ‘quantreg’package (Koenker, 2015) available in R software (R Develomment Core Team, 2015). This robust method is especially useful for estimating rates of change along or near the upper limit of the functional size-density relationship (Cade and Noon, 2003; Zhang et al., 2005). Charru et al. (2012) demonstrated the suitability of National Forest Inventory data for estimating self-thinning relationships, based on the assumption that a fraction of the plots are located in stands undergoing density-dependent mortality over the observed species mean diameter range. The logarithmic form of equation (3) was used in this analysis, as follows: ββ ε=+ + ()Ndgln ln 4 01 where β 0 and β 1 are parameters of the self-thinning line to be estimated. We focused attention on upper quantiles 90–99 th .Theparameters obtained in the outermost quantiles were compared with the slope and intercept fitted by ordinary least squares regression (OLS). Additionally, we tested whether other independent variables such as stand origin might significantly affect the intercept value in the self-thinning line. Mixing effects at stand level The volume increment per hectare (iV) between the two inventories was divided by the time span between them, because the period between measurements varied among plots (10–13 years). The growth efficiency = I Vp iV P/ ii i , was defined as the volume increment of species iV i (m 3 ha −1 year −1 ) divided by the stocking proportion (equation 1) or area occupied by the species. To calculate total stand density ( SD ) for comparing growth in pure and mixed stands, the relative total stand density 100 100 Kilometers50 0 N E S W Figure 1 Location of selected NFI plots and the distribution of P. pinaster (gray) and P. sylvestris (dark gray) in Spain. Mixing effects on growth efficiency in mixed pine forests 383 Downloaded from https://academic.oup.com/forestry/article/90/3/381/2452295 by UNIVERSIDAD DE VALLADOLID user on 23 November 2021 was used based on the same concept as the stocking proportion (Sterba et al., 2014): =+ ()SD SDIR SDIR 5 ij Similarly, stand variables such as H oand d gwere included in the model as surrogate variables to account for the high variability in the site qualities and ages presented in NFI data. To test if the species mixture affects growth efficiency and study the density-growth relationships, the following general model was fitted for Scots pine and maritime pine, respectively: ε () =+ ()+ () +()++ + () IVpaaHoadgaSDaPaln ln ln ln ORIG 6 ii ji01 2 3 4 5 where the dependent variable is the volume growth efficiency of species i. Dominant height (H O ) of the stand, the quadratic mean diameter ( d g) of the species i and total stand density ( SD ) are the independent variables. To test if species mixture affects growth, we included the species proportion by area of the complementary species j ( P j ). A positive coefficient for this variable would indicate that growth efficiency of the species i increases with an increasing proportion of the species j , positive mixing effect on the growth of the analyzed species. Additionally, a dummy variable for each origin was included ORIG (0 for natural plots or 1 for plantation plots) and interactions between species proportion by area and independent variables were tested. Finally, we tested the benefit of including the set of climate variables in the model structure. Linear regression was fitted using the least squares approach establishing significance at P<0.05. The relevance of the species proportion in the growth efficiency model was tested using the F-statistic to compare equation (6) (mixture model) against a null model (without the P j parameter). Mixing effects on inter-tree competition Single-tree models are a very effective tool for exploring the effect of competition on growth. Despite the uncertainty of using expansion factor to estimating stand variables in NFI plots (trees are recorded according to dbh class over circular concentric plots), these coefficients are able to produce unbiased values (Bravo et al., 2002). Thus, a distance independent approach is required to quantify competition with this kind of data. Stand variables such as basal area, stand density, number of neighbouring trees, basal area of trees larger than the target tree or size of individual trees relative to mean tree size (Larocque et al., 2013) can be used for this purpose. With the same plots used to assess mixture effects at the stand level, a tree-level dataset was selected that incorporated the presence of both species in the plot regardless of their mixture proportions. Dead trees in either inventories or plots with anomalous growth data were Table 1 Main stand variables for mixed and pure plots at the beginning of the growth period (used to fit the stand-growth models) Origin Mixed plots Pure plots Total stand Scots pine Maritime pine Scots pine Maritime pine NPNPNPNPNP No of plots 82 33 82 33 82 33 165 45 106 42 N(trees ha −1 ) Mean 564.5 1023.6 309 590 264 423 667.9 1293.3 796.5 863.6 Max 1479.9 3087.5 1322.7 1987.6 1177.7 1481.9 2893.0 3755.9 2291.8 1955.8 Min 87.1 226.4 14.2 88.4 19.2 28.3 65.1 240.5 104.1 283.0 SD 312.5 744.6 241.4 500.1 218.1 398.5 468.6 678.6 558.3 450.7 G(m 2 ha −1 ) Mean 26.8 21.6 12.7 9.7 14.3 11.7 31.9 21.2 30.3 23.7 Max 53.7 41.4 35.4 22.9 44.8 31.9 68.9 54.0 67.4 48.0 Min 8.4 7.9 1.5 1.5 0.8 1.9 8.9 7.4 8.0 8.1 SD 12.3 9.5 8.3 5.3 9.3 7.5 14.2 10.7 15.5 10.3 H o (m) Mean 14.0 9.8 13.3 8.9 13.3 9.8 15.4 8.7 12.0 9.7 Max 22.3 13.6 21.0 14.1 21.3 13.2 27.9 19.2 20.1 14.3 Min 6.3 5.3 5.8 4.2 4.8 5.6 6.1 4.7 5.1 5.5 SD 3.4 2.2 3.5 2.6 3.8 2.0 4.0 3.2 3.2 2.3 dg (cm) Mean 26.3 18.0 25.1 17.0 28.9 21.5 27.7 15.1 24.3 19.6 Max 39.5 33.6 50.7 35.8 41.6 31.7 51.5 32.5 46.1 33.5 Min 11.7 10.7 9.8 8.8 10.3 11.2 10.5 8.9 10.0 10.7 SD 7.0 5.3 7.8 6.8 8.9 5.5 8.7 4.6 7.6 5.0 Volume increment (m 3 ha −1 yr −1 ) Mean 5.1 7.4 2.3 3.2 3.1 4.1 4.6 6.7 5.0 7.7 Max 12.5 17.6 8.5 7.5 9.2 14.4 14.0 16.6 14.4 19.1 Min 0.3 1.4 0.1 0.6 0.1 0.6 0.1 1.6 0.7 1.0 SD 2.8 4.1 1.7 1.8 2.3 3.2 3.1 3.5 3.1 4.8 Total stand density (SD) Mean 0.48 0.41 0.23 0.19 0.25 0.22 0.55 0.41 0.56 0.44 Max 0.99 0.76 0.62 0.47 0.81 0.60 1.17 1.00 1.2 0.88 Min 0.17 0.16 0.03 0.03 0.01 0.03 0.15 0.15 0.15 0.15 SD 0.23 0.18 0.15 0.11 0.17 0.14 0.24 0.19 0.28 0.19 N, stems per hectare; G, basal area per hectare; H o , dominant height; dg, quadratic mean diameter; N, Natural stands; P, Plantation stands; SD, Total stand density as the sum of the relative stand density of each species SD =SDIR i +SDIR j . Forestry 384 Downloaded from https://academic.oup.com/forestry/article/90/3/381/2452295 by UNIVERSIDAD DE VALLADOLID user on 23 November 2021 excluded. Individual tree data for both species at the beginning of the studied growth period is shown in Table 2. An adequate individual tree growth model must include both an expansion and decline component to represent the growth pattern (Zeide, 1993). Equation (7) includes both components as a function of tree size. The competition status of the tree was included as a modifier component ( e Comp ) that reduces the potential growth rate to the actual growth rate (Bravo et al., 2001;Pretzsch and Biber, 2010;Río et al., 2014): =⋅ ⋅ ⋅⋅ ()iv a hedbh e 7 iai a i a 0dbh Comp i 123 where iv i is the volume increment of a tree with a given diameter at breast height dbh i , and total height h i ; a 0 ,a1,a2and a 3 are species-specific parameters and Comp is the competition term. We used the stand density index concept ( S D I in equation 2) as a measure of size-symmetric competition (S). To analyze size-asymmetric competition (AL), we used the stand density index of trees larger than the target tree ( S DIL ) (equation 8) as an indicator of asymmetric competition for light (Pretzsch and Biber, 2010): ⎛ ⎝ ⎜⎞ ⎠ ⎟ =()Ndg SDIL 25 8 L L E where NLand d g L , are the number of stems per hectare and quadratic mean diameter of trees larger than the target tree, respectively, and E is Reineke’s maximum stand density exponent. We used species-specific values of E for each species, which were fitted using the QR technique (see ‘Stocking proportion and maximum stand density’section and Table 3) instead of the generic value proposed by Reineke (1933). Both size-symmetric and size-asymmetric competition were calculated in relative terms (equation 2) to take into account species-specific site occupancy, as the ratio between S D I and S DIL , respectively, and the maximum stand density index of each species (Sterba et al., 2014; Condés and Río, 2015). To distinguish between intraspecific and interspecific competition and evaluate the influence of the mixture on tree growth, the relative size-symmetric ( S DIR ) and relative sizeasymmetric (SDIRL) competition indices attributed for each species were calculated, following a similar approach to that of Río et al. (2014). Resulting in six models for each species, three size-symmetric/sizeasymmetric structures (S, AL, S +AL) by the combination of two species composition structures (total, intra +inter). If inclusion of the total stand density indices in the growth model results in a better fit than the use of the intraspecific and interspecific indices, it may indicate similar intraspecific and interspecific competition or no mixing effect (Río et al., 2014). The performance of the different competition structures was analyzed in the tree volume growth model (equation 7), using informationtheoretic approach (Anderson, 2008), which allows evaluate multiple non-nested models relative to each other and quantify the relative support for multiple models simultaneously (Zuur et al., 2009), especially useful when effects of different competition structures on tree growth are assessed (Papaik and Canham, 2006). Fitted models were ranked by their AICc value (Second-order Akaike Information Criterion), model with the lowest AICc could be considered most parsimonious ‘best model’. The absolute magnitude of the differences in AICc between alternate models (ΔAICc) and Akaike weights (ωi) provides an objective index of the strength of empirical support for the competing models or alternate hypothesis (Anderson, 2008). The multi-model inference procedure was performed with MiMIn R-package (Bartó, 2016). To attain normal distribution of the residuals and reduce heteroscedasticity, a natural logarithmic transformation of the original dependent variable was applied. Given the hierarchical nature of the data, we used a linear mixed-model approach (equation 9). We tested random effects at two nested levels: phytoclimatic region and plot-in-region. However, the inclusion of the phytoclimatic region in the hierarchical structure did not improve the models performance; thus, we only used plot as the grouping structure of the random effects. We included random effects in both intercept and logarithm diameter terms. The inclusion of random effects allows relaxing the assumption of independence of residuals as Table 2 Tree characteristics and competition status by stand origin in selected plots at the beginning of the growth period (used to fit tree-growth models) Origin ndbh (cm) h(m) iv (dm 3 yr −1 ) SDIR SDIRL SDIR intra SDIR inter SDIRL intra SDIRL inter Pinus sylvestris N 1317 Mean 27.2 14.1 11.12 0.52 0.30 0.33 0.19 0.16 0.14 Min 7.5 3.0 0.01 0.09 0.00 0.01 0.01 0.00 0.00 Max 54.1 24.0 68.15 1.06 1.04 0.97 0.94 0.96 0.94 SD 9.12 3.93 8.31 0.23 0.21 0.21 0.17 0.15 0.15 P 419 Mean 17.5 8.9 9.16 0.36 0.23 0.19 0.18 0.09 0.14 Min 7.5 3.5 0.20 0.07 0.00 0.01 0.01 0.00 0.00 Max 55.8 18.0 74.89 0.65 0.62 0.39 0.49 0.35 0.49 SD 6.71 2.48 9.23 0.16 0.14 0.10 0.11 0.08 0.10 Pinus pinaster N 1367 Mean 32.3 14.0 18.69 0.54 0.23 0.38 0.16 0.18 0.05 Min 7.65 2.0 0.01 0.09 0.00 0.01 0.01 0.00 0.00 Max 73.2 25.5 120.51 1.06 0.98 0.94 0.97 0.92 0.44 SD 11.4 4.40 14.50 0.25 0.19 0.24 0.14 0.17 0.07 P 523 Mean 22.3 9.6 17.53 0.39 0.15 0.26 0.13 0.13 0.02 Min 7.7 4.0 0.15 0.07 0.00 0.01 0.01 0.00 0.00 Max 54.5 19.0 76.54 0.65 0.63 0.49 0.39 0.48 0.26 SD 6.91 2.10 11.49 0.16 0.12 0.13 0.09 0.11 0.04 dbh, diameter at breast height; iv, tree volume growth; SDIR, size-symmetric competition index; SDIRL, size-asymmetric competition index for larger trees; intra, intraspecific competition; inter, interspecific competition; SD, standard deviation; min, minimum value; max, maximum value; N, Natural stands; P, Plantation stands. Mixing effects on growth efficiency in mixed pine forests 385 Downloaded from https://academic.oup.com/forestry/article/90/3/381/2452295 by UNIVERSIDAD DE VALLADOLID user on 23 November 2021 they would absorb the unaccounted variability at the aforementioned levels due to growing conditions other than competition: ε ( )=( + )+ ⋅ +( + )⋅ ( ) +⋅()+ ∑++ () iv a u a a u ah b a ln dbh ln dbh ln Comp ORIG 9 ij j i j i im im ij 00 1 22 34 where iv ij is the volume increment of the i th tree measured in the j th plot, u j0and u j2are plot level random effects and ε ij the error term. Comp im represents the m different terms of the competition for tree i according to the different structures mentioned and b m their corresponding parameters to be estimated. ORIG is a dummy variable for the origin of the stands. All errors terms are assumed to be normally distributed with mean 0, σσεσ( )~( ) ( )~( ) ~( ) uN uN N0, , 0, and 0, jj ij00 222 2res 2. To determine the effect of the nested random structure in the model we used the Akaike’s information criterion (AIC) and log-likelihood estimates for comparison. We assumed that the random effects are independent; in this case we used a diagonal form for variance–covariance matrix. Models were fitted using functions from the nlme R-package (Pinheiro and Bates, 2000) and the maximum likelihood method to allow comparison among different competition structures. The unbiased restricted maximum likelihood (REML) method was then used to fit the parameter estimates and variance components for the final model that was selected for each tree species. For each model, we calculated conditional R 2 values, which account for the explanatory power of both fixed and random effects and marginal R 2 to describe the proportion of variance explained by the fixed factors alone, both as a measure of goodness-of-fit of the models (Nakagawa and Schielzeth, 2013). Results Self-thinning boundary line and maximum stand density A maximum size-density relationship was clearly identified based on the NFI data (Figure 2). No effects from stand origin were detected in the self-thinning limit parameters, so a unique model was fitted for each species. We chose the 0.95 quantile because it showed significant parameters and the closest 95% confidence intervals of the highest quantiles. Additionally, the parameters of both species in this quantile were different than the OLS estimation, which represents the central tendency line across the range of data (Supplementary data, Figure S2). The intercept and slope of the self-thinning relationship for both species are given in Table 3; these parameters were used in equation (2) to estimate the stand density and S DI max values. Growth efficiency and effect of species mixture at stand level The estimated parameters and regression statistics of the growth efficiency models that were found for the two species using the Table 3 Quantile regression parameters and standard errors (in parentheses) for linear models of self-thinning relationships of pure stands (Eq. 4) Species Intercept Slope SDI max Scots pine 12.801 (0.192) −1.789 (0.061) 1143.65 Maritime pine 13.218 (0.278) −1.929 (0.088) 1103.64 SDI max , maximum stand density index estimated using dg =25. Pinus sylvestris Pinus pinaster 2000 4000 6000 2000 4000 6000 20 30 40 dg (cm) N (trees ha–1) Figure 2 Self-thinning boundary line using NFI data plotted on a log–log scale. Lines by lines joining dots represent the ongoing density-related mortality within plots between the second and third NFI. Straight lines represent maximum size-density lines (solid =QR, dotted =function fitted by Río et al., 2001,2006). Forestry 386 Downloaded from https://academic.oup.com/forestry/article/90/3/381/2452295 by UNIVERSIDAD DE VALLADOLID user on 23 November 2021 stocking proportion approach are shown in Table 4. All coefficients were significant in the final models (P<0.05) and a similar relationship between growth efficiency and predicted variables was found in both species: a positive relationship with stand dominant height, total stand density and stocking proportion of complementary species, but a negative relationship with the quadratic mean diameter of the target species. Dominant height and quadratic mean diameter had been used as surrogate variables of site quality and age, respectively; thus, the results in Table 4mightindicateanincreasein growth efficiency with increasing site quality and a decrease of growth efficiency with aging. Growth efficiency increased with total stand density and obtained maximum volume increment in fully stocked stands (SD =1), as expected. The species mixture had a positive and significant effect (P<0.05) on the growth efficiency of both species (Table 4) and significant difference with the null model revealed the influence of species proportion in the mixture model, indicating that both species grew better in mixed stands than in pure stands. The interaction between stand density and stocking proportion of species by area tested did not improve the growth efficiency models in either species. Additionally, the stand origin was statistically significant, indicating that a model is required for each origin and species. Inclusion of climatic variables did not prove to be statistically significant. Comparison of tree-level competition structures in mixed stands Table 5compares the different models that include sizesymmetric or size-asymmetric competition structures (S, AL, S + AL) and assumes different composition in the competition structures (both species together, intraspecific and interspecific). Models that simultaneously included size-symmetric and sizeasymmetric terms (S +AL) outperformed the models than used only S and AL in both species composition structures. Further, models expressed only by size-asymmetric competition (AL) resulted better and the parameter was larger than size-symmetric competition (S), which denotes the importance of light competition in these species. The combination of size-symmetric and size-asymmetric terms (S +AL) with different species composition (intra+inter) emerged as the best competition structure model, with lowest AIC c values. Contrasting results were observed between species when best competition structures were compared. In the case of Scots pine, intraspecific and interspecific competition were significant in size-asymmetric competition terms with higher coefficient value of intraspecific term that means that interspecific competition was less intense for Scots pine. In size-symmetric competition the variable related to maritime pine competition was found to be non-significant, which means no effect of interspecific competition on growth. For maritime pine, the same pattern was observed in sizesymmetric competition structure than in Scots pine; instead in size-asymmetric competition the coefficient related with intraspecific parameter was negative and lower than the Scots pine competition term. However, this model was slightly better (ΔAIC c =0.44) than the second best ranked model (without different species composition in the competition). Despite the uncertainty to define the ‘best model’, both models indicate no benefit from mixture at least by size-asymmetric competition, which is more intense for tree growth than that of sizesymmetric competition. The volume growth models with lowest AIC c value were fitted with the REML procedure (Table 6). Parameters were significant at 95% probability. Figure 3depicts the effects of size-symmetric and size-asymmetric competition analyzed independently, when the target tree was influenced by intraspecific competition only (pure stand) and when competition was the combination of intraspecific and interspecific interactions (mixed stand). In order to illustrate the different values of competition structures, simulations of tree volume growth were set on an average tree with 30 cm dbh and 14 m hin a natural stand and using mean competition status values (Table 2). The stronger the size-symmetric interspecific competition, the higher the profit from the mixture on both species (Figure 3,above). Positive mixture effect on tree growth was evident for Scots pine (Figure 3, below), while the slight difference with model without splitting the competition by species highlighted the strongly effect of intraspecific size-asymmetric competition in maritime pine. Discussion Feasibility of size-density relationship estimation with NFI data Appropriate estimation of potential or maximum density is a key step in defining the stocking proportion of species in mixed stands. The attributed mixing effects introduced by the use of an inadequate species proportion definition can cause confusions in Table 4 Estimated coefficients and P-values (P) for volume growth efficiency models (equation 6) of the two species (standard errors in parentheses) Scots pine Maritime pine Coeff. PCoeff. P Intercept 3.516 (0.461) <0.001 3.782 (0.475) <0.001 ln (H o ) 0.626 (0.244) 0.010 0.530 (0.209) 0.011 ln (dg i )−1.076 (0.176) <0.001 −0.971 (0.150) <0.001 ln (SD) 0.238 (0.117) 0.042 0.331 (0.090) <0.001 P PT 0.341 (0.148) 0.021 P PS 0.355 (0.149) 0.018 Origin i −0.306 (0.122) 0.013 −0.338 (0.094) <0.001 MSE 0.541 0.389 RSE 0.735 0.623 Adj.R 2 0.222 0.228 F5.12 0.024 5.63 0.018 ∆ AIC 3.2 3.7 H o , Dominant height; dg, quadratic mean diameter; SD, total stand density; P PT , stocking proportion of Maritime pine; P PS , stocking proportion of Scots pine; Origin, dummy variable (0 =natural stands, 1 =plantation stands); RSE, residual standard error; Adj.R 2 , adjusted coefficient of determination R 2 ;F,F-statistic and ΔAIC, delta AIC between mixture versus null model. Mixing effects on growth efficiency in mixed pine forests 387 Downloaded from https://academic.oup.com/forestry/article/90/3/381/2452295 by UNIVERSIDAD DE VALLADOLID user on 23 November 2021 the interpretation of growth comparisons between mixed and pure stands (Sterba et al., 2014). QR models reached lower maximum size-density line and S DI max value than the upper limiting boundary lines fitted in previous studies for these species in the region (Río et al., 2001, 2006) (Figure 2), mainly due to differences in the intercept values (C), which may vary with site conditions. This was reported in other studies where larger intercept values were found in stands growing on more productive land (Bi, 2004; Weiskittel et al., 2009), causing considerable differences in density. Though stand origin or purity (Weiskittel et al., 2009)have also been associated with significant changes in the intercept on the self-thinning boundary line; however, we did not observe effect of stand origin in our results. The similarity of species-specific slope observed for Scots pine (Supplementary data, Figure S2) with the parameters proposed by Río et al. (2001) estimated from non-thinned permanent plots, suggesting that the QR approach is appropriate for estimating the self-thinning boundary line using inventory data. In the other hand, the self-thinning coefficient found for maritime pine in this study support the idea of using a species-specific slope (Charru et al., 2012) rather than a constant coefficient among species. However, Río et al. (2006) used the slope stated by Reineke to fit the self-thinning line for P. pinaster, which did not allow us to compare and determine whether our selfthinning relationship was correct. Effects of the species mixture on the growth of Scots pine and maritime pine stands Site quality is a trigger factor and could determine the effect of the mixture (Forrester, 2014;Huber et al., 2014). Though stand age is missing in NFI data and consequently site index could not be estimated in this study, we used dominant height and quadratic mean diameter as surrogate variables for site quality and stand age or development stage, respectively, in the growth efficiency model. The traditional approach to calculate side index based in Eichorn’s rule might not fit in mixed stands due to interspecific interactions could modify the expected side index from pure stands (Río et al., 2016). However, site index for pure stands (Pretzsch et al., 2015b) or abiotic gradients to describe site productivity (Toïgo et al., 2015) are commonly used to analyze the relationship between total yield and site conditions in mixed forests. In this study, we understand that at a given Table 5 Ranking of growth models comparing competition structures and their parameters estimated. The competition status was defined as the relative stand density index calculated for size-symmetric (SDIR) and/or size-asymmetric competitors (SDRIL) and for trees of all species or divided by intraspecific and interspecific competitors SDIR SDIRL SDIR inter SDIR intra SDIRL inter SDIRL intra AIC c ΔAIC c ω i Pinus sylvestris ns −0.8558 −0.8312 −1.0079 3077.7 0 0.801 −0.4552 −0.8943 3081.0 3.25 0.158 −1.0860 3083.6 5.91 0.042 −0.8985 −1.1969 3104.6 26.92 0.000 −0.7641 −1.3640 3118.8 41.08 0.000 −0.9990 3120.7 43.04 0.000 Pinus pinaster ns −0.6262 −1.614 −0.7307 3202.9 0 0.535 −0.5251 −0.8650 3203.3 0.44 0.429 −1.669 −0.9077 3209.0 6.12 0.025 −1.0120 3210.6 7.68 0.012 −0.8950 3233.1 30.17 0.000 −0.7357 −0.9379 3234.9 32.02 0.000 SDIR, size-symmetric competition index; SDIRL, size-asymmetric competition index for larger trees; intra, intraspecific competition; inter, interspecific competition; AICc, Second-order Akaike Information Criterion; ΔAIC c , difference between the best model and the i th model; ω i , Akaike weights. Significance of coefficients P≤0.05; ns, not significant. Table 6 Estimated coefficients for the best tree growth models (equation 9) of the two species (standard errors in parentheses), including competition structures. Conditional R 2 values are showed as a measure of goodness-of-fit Scots pine Maritime pine Intercept −4.0212 (0.4170) −3.0889 (0.3921) dbh −0.0465 (0.0091) −0.0242 (0.0062) ln (dbh) 2.2426 (0.2089) 1.8626 (0.1780) ln (h) 0.3951 (0.1050) 0.3854 (0.0935) Origin −0.3493 (0.0894) −0.3442 (0.0821) SDIR intra −0.8606 (0.2638) −0.6377 (0.2089) SDIR inter −0.2275 ( 0.314) ns −0.401 (0.2509) ns SDIRL intra −1.0020 (0.2053) −0.7242 (0.1629) SDIRL inter −0.8196 (0.3204) −1.6234 (0.3853) AIC 3104.13 3230.74 R 2 conditional 0.62 0.60 R 2 marginal 0.43 0.46 Variance of random effects Intercept 0.3637 0.3229 ln (dbh) i 8.479 ×10 −7 0.0403 Error 0.0003 0.0086 R 2 conditional, variance explained by both fixed and random effects; R 2 marginal, describes the proportion of variance explained by the fixed factors alone. Significance of coefficients P≤0.05, ns, not significant. Forestry 388 Downloaded from https://academic.oup.com/forestry/article/90/3/381/2452295 by UNIVERSIDAD DE VALLADOLID user on 23 November 2021 quadratic mean diameter, dominant height can be a descriptor for site quality, with the constrain that dominant height cannot be directly associated to site quality class, so, mixture effect at stand level throughout a productivity gradient could not be assessed (Condés et al., 2013). The criteria used to select the pure and mixed plots allowed us to reduce the shift in ecological variables—between pure and mixed stands—and the effects associated with silvicultural practice that control the level of growing stock that might modify (Forrester and Pretzsch, 2015) or mask the real mixture effect. However, the low variability explained by the growth efficiency models in both species may be due to the uncertainty related with NFI data characteristics and methodological approach: stand age information is unavailable, uncertainty about past silvicultural treatments, and estimated stand variables using expansion factors could generate unbiased but uncertain errors (Bravo et al., 2002). Further, the self-thinning law is difficult to verify and we assume that stand density in mixed stand can be estimated accordingly by relative densities of both species, this implies risk to include bias in the estimated yielding in mixed stands. In this sense, when highly accurate predictions are required, some methods exist to deal with bias in predictions, such as marginalization of predictions over the error term distribution through numerical integration (Fortin, 2013). Nevertheless, the method to estimate species occupancy by area in mixed stands and therefore species proportion has been compared with other methods and has provided often good results (Huber et al., 2014;Sterba et al., 2014). Spanish NFI data offered the advantageous possibility of obtaining a general overview of the behaviour of this mixture type in an unbiased and systematic sample distributed across the complete range of mixture proportions. Thus, we compared the ‘true mixture effect’(Forrester and Pretzsch, 2015) independently of the mixture proportion, revealing apparent profitatthestand scale and mutual benefit of both species from the mixture. Given that complementarity processes require niche differentiation among species (Kelty, 2006), mutual facilitation in which both species gain (Larocque et al.,2013), as we observed in our results, is an unexpected indication of a complementary relationship between Scots pine and maritime pine. When light demanding species with similar crown architectures coexist such as these two pines, complementarity is less likely to occur. Our results suggest that even small differences in species traits, as small differences in shade tolerance, can trigger positive interspecific effects, through mechanisms of competition reduction or complementarity. Species mixture is able to induce changes in crown dominance, above-ground species allometries and structural stands traits (Pretzsch and Schütze, 2005;Pretzsch et al.,2016)which are related with the way light is distributed among trees and hence productivity. This has sensitive effects for light demanding species mixture (Jucker et al.,2014). The potential growth gain of Scots pine was also found in mixture with Quercus pyrenaica (Río and Sterba, 2009) and Fagus sylvatica (Condés et al., 2013;Pretzsch et al., 2015b), and Pinus nigra (Jucker et al., 2014), but no such effect was reported in mixture with Quercus petraea (Toïgo et al., 2015). Most of these cases, complementary strategies of light interception or light use efficiency between species were the main driver for the positive mixing effects. To our knowledge, no prior research has reported mixture effects on growth efficiency for maritime pine. Despite mixture effect observed in this study, the mechanisms that trigger species interactions are not easy to disentangle in a stand-level analysis, because some factors may act simultaneously. In climatically stressful Mediterranean conditions, there is evidence of higher wood production in mixed forests due to species niche partitioning in water use (Vilà et al., 2007). Competition effects on tree growth Size-symmetric and size-asymmetric competition have frequently been related to below-ground and above-ground resources, respectively (Larocque et al., 2013). Differentiation of intraspecific and interspecific competition has also provided useful information for inferring possible factors influencing tree growth (Manso et al., 2015) and mortality (Condés and Río, 0 5 10 15 20 25 30 0 0.2 0.4 0.6 0.8 1 iv (dm3 yr –1)iv (dm3 yr –1) Symmetric competition 0 5 10 15 20 25 30 0 0.2 0.4 0.6 0.8 1 Asymmetric competition PS intra PS intra+inter PT intra PT intra +inter Figure 3 Volume growth as a function of size-symmetric competition (top, asymmetric competition =0) and size-asymmetric competition (bottom, symmetric competition =1) in Scots pine (PS, circles) and maritime pine (PT, triangles) for a tree with 30 cm dbh and 14 m hin a natural stand and using mean competition status values (Table 2). Solid lines indicate results for intraspecific competition only (pure stand); dashed lines indicate results for a combination of intraspecific and interspecific competition structures (mixed stand). Mixing effects on growth efficiency in mixed pine forests 389 Downloaded from https://academic.oup.com/forestry/article/90/3/381/2452295 by UNIVERSIDAD DE VALLADOLID user on 23 November 2021