Maximum stand density strongly depends on species-specific wood stability, shade and drought tolerance

Maximum stand density strongly depends on species-specific wood stability, shade and drought... Abstract The abundance of stems in crowded populations and the subsequent self-thinning is a key issue in forest stand dynamics. However, the mechanisms that control self-thinning are challenging to model. Although some attempts to include climate and structural traits like specific gravity (SG) are promising, they remain confined to North American species. In this study we aimed to disentangle how SG along with two major abiotic stress tolerances, i.e. shade and drought tolerance, contribute to the maximum density of a forest stand across a climatic gradient in Europe, and thus test the validity of the species-specific trait control over stand density. We propose a modelling approach that incorporates the tolerance to drought and shade in the determination of maximum relative stocking. Here, relative stocking refers to the degree of tree crowding in forest ecosystems. A relative stocking base model where specific wood density is inversely related to stand density is modified, adding normalized indices of drought and shade stress tolerance. We used available species tolerance rankings modulated by stress intensity to analyse the effects of abiotic stress polytolerance or trade-offs in the study area which represent an environmental gradient from Alpine to Mediterranean climate in northern Spain. Results indicated that the role of drought tolerance in controlling maximum stand density is stronger in warmer sites. The simultaneous tolerance to shade and drought results in less carrying capacity of sites. In those sites where there is no water limitation but minimum temperature is very low the tolerance to bending stress (i.e. specific gravity) explains better the maximum tree occupancy. Introduction Forest stand density is a key variable in forestry and ecology. Understanding the reasons why tree species occupy an area and the maximum capacity of each site to harbour a certain number of individuals is at the core of management decisions. Early silvicultural studies highlighted the concept that stand density and its control allows managers to keep production within optimum levels to compensate the growth trade-off between stand- and tree-level (Langsaeter, 1941; Assmann, 1971). Long (1985) summarized this hypothesis establishing the density limits between the thresholds of lower full site occupancy and the onset of self-thinning. These limits are based on the maximum stand density of a site that represents the maximum degree of site occupancy (Dean and Baldwin, 1996). The determination of maximum density has largely been explored in terms of maximum size-density relationships (MSDR). Most forestry studies explore MSDR using Reineke’s equation (Reineke, 1933) and the associated Stand Density Index (SDI) that relates the number of stems as a function of the size of the trees in relation to a fixed diameter, usually 25 cm (Pretzsch and Biber, 2005). Similarly, ecological studies have dealt with the −3/2 self-thinning law postulated by Yoda et al. (1963) that relates the individual mean volume with the number of individuals. Both relationships in their original formulation are based on the assumption that there is a limit in the relationship between size and density and they are equivalent if tree volume is proportional to quadratic mean diameter (Burkhart and Tomé, 2012). Maximum forest stand density represents an upper limit to the occupancy of a site, and growth is only possible after the death of some individuals. This upper limit on potential site occupancy has been considered to be species- and site-specific by several authors using different variables to characterize the stand. Condés et al. (2017) used de Martonne’s aridity index (de Martonne, 1926) for Scots pine (Pinus sylvestris L.) and European beech (Fagus sylvatica L.) whereas Zhang et al. (2013) included site index to model MSDR for ponderosa pine (Pinus ponderosa). The models are fitted for specific species compositions and the results cannot be applied to different compositions. Thus, while potentially useful for well-defined management situations, the ecological inferences from such models remain rather limited. The simultaneous representation of density boundary lines across many different species, growing in both pure and mixed stands and in different environments is still seldom explored. Site index was also used in conjunction with other stand-specific attributes like site index, climate, stand purity and origin to fit a MSDR for different forest populations and communities of Douglas-fir (Pseudotsuga menziesii var. menziesii [Mirb.] Franco), western hemlock (Tsuga heterophylla (Raf.) Sarg.) and red alder (Alnus rubra Bong.) (Weiskittel et al. 2009). Similarly, site quality and species composition were included in modelling maximum density of trembling aspen (Populus tremuloides Michx.) and white spruce (Picea glauca (Moench) Voss.) stands (Reyes-Hernandez et al., 2013). Another approach is that based on species-specific traits like specific wood gravity, SG (Woodall et al., 2005; Ducey and Knapp, 2010) or a combination of two traits, like SG and shade tolerance (ST) modulated by climate (Ducey et al., 2017). The hypothesis behind the use of species-specific traits to predict maximum density is based on stress tolerance. Specific wood density is associated with tolerance to bending stress, as developed by Dean and Baldwin (1996) who first postulated that species with light wood would need more stems to support the same foliage mass per unit area than species with harder wood. The inverse relationship between maximum stand density, defined as Reineke’s stand density and SG was presented by Woodall et al. (2005) whose relationship can be applied to mixtures as long as specific gravity is known for all component species. Ducey and Knapp (2010) combined Woodall et al.’s (2005) approach with the summation method to apply the additive stand density index (ASDI) developed by Long and Daniel (1990) in heterogeneous forests along with the Curtis’ (1971) modified version of the tree area ratio (Chisman and Schumacher, 1940). The resulting model predicts the relative stocking of a stand based on the specific gravity of tree species. Recently, Ducey et al. (2017) added tolerance to low light environments in addition to bending stress. Their results showed how climate variables modulated stress tolerance and provided the improvement to a more mechanistic explanation of site occupancy. Shade tolerance is a key mechanism shaping forest communities (Valladares and Niinemets, 2008) and it has been used to estimate maximum stand density (Ducey et al., 2017). Other abiotic stresses might be more important, however, in other climatic regions like the Mediterranean where drought events are more harmful than low light levels. Thus, we must test the idea that a single or two-trait approach may be insufficient to characterize maximum density and raise the issue of dealing with polytolerance, i.e. the ability of a species to simultaneously tolerate different abiotic stressors. In this paper we aimed to identify which trait contributes more to the maximum density of a forest stand, and if the relative contribution of two major types of abiotic stress tolerance, i.e. shade and drought tolerance, changes across a climatic gradient. Our working hypotheses are: Mechanical properties, like the wood specific gravity, have an inverse relationship with stocking of forest stands across a climatic gradient. Shade tolerance value, as modulated by cold temperature, affects maximum density inversely: the more shade tolerant the species and colder the site the less number of individuals needed to occupy the stand. Drought tolerance is more important in drought-prone areas, like in Mediterranean areas, where fewer individuals are needed to fully occupy the stand due to site conditions. Methods Study area, stand and climate data The study is carried out in a transition zone covering Mediterranean, Atlantic, and Alpine climate in Navarra province, northern Spain (Figure 1). The mean annual temperature is 10.5°C, ranging from 4.3 to 14.2°C, and annual rainfall ranges from 373 to 2629 mm, averaging 1378 mm (Gonzalo-Jiménez, 2010). Stand information comprising species composition, number of individuals by unit area (N, stems ha−1), and quadratic mean diameter (dg, cm) were extracted from the Third Spanish National Forest Inventory (NFI) available at www.mapama.org. Each NFI plot consists of four concentric subplots with radius 5, 10, 15 and 25 m where trees with DBH equal to or over 7.5, 12.5, 22.5 and 42.5 cm were measured, respectively. To get the number of stems by unit area an expansion factor is applied for each subplot, i.e. for the 10 m subplot each tree represents 10 000/(π × 102) individuals. We removed from the database plots in fertilized or irrigated Populus x canadensis plantations to make sure that intensive management did not affect the intensity of abiotic stress. We retained 2803 plots in four climatic regions built upon Allué-Andrade’s (1990) phytoclimatic classification: Mediterranean, Mountainous, Maritime and Transition (see Supplementary Table 1 for description of climatic regions). A species-specific wood gravity for each species were retrieved from Gutiérrez and Plaza (1967) whereas shade and drought tolerance (ST and DT, respectively) was obtained from Niinemets and Valladares (2006). For those species without a tolerance value we applied the average value for the genus. Table 1 shows main stand level properties and stress tolerance range by climatic gradient. The dataset comprises 66 species ranging from temperate to Mediterranean-type forests. The most abundant species are European beech (Fagus sylvatica L), Scots pine (Pinus sylvestris L.), Holm oak (Quercus ilex L.), oak (Q. robur and Q. petraea) and quejigo oak (Quercus faginea Lam.) (Supplementary Table 2 describes the main dendrometric characteristic of each species growing in the study region). Figure 1 View largeDownload slide Study area (Navarra, Spain), biogeoclimatic zones and National Forest Inventory plots used in this study. Figure 1 View largeDownload slide Study area (Navarra, Spain), biogeoclimatic zones and National Forest Inventory plots used in this study. Table 1 Descriptive stand variables, specific gravity and stress tolerances by bioregions in the study area. Bioregion # plots Hdom (m) BA (m2 ha−1) dg (cm) N (stems ha−1) SG ST DT Maritime 1044 17 (5.3) 3.5–33 28.1 (10.6) 0.7–78.2 26.1 (11.5) 8–111.5 755.6 (577) 5.1–3738.4 0.7 (0.1) 0.4–1 3.4 (1.1) 1.2–4.6 3 (0.8) 1.2–5 Mediterranean 374 12.4 (6) 3–27.6 22.8 (11.5) 0.4–52.6 19 (7) 7.6–81.3 998.8 (677.9) 10.2–2910.8 0.7 (0.2) 0.4–1 2.4 (0.8) 1.4–4.6 3.6 (1.1) 1.2–5 Mountainous 1121 17.8 (6) 2.5–36.5 31.5 (13.8) 0.4–82.7 25.6 (10) 7.6–95.5 808.2 (561.8) 5.1–3243.2 0.7 (0.1) 0.4–1 3.2 (1.4) 1.2–4.6 3.3 (0.9) 1.7–5 Transition 264 10 (4.1) 3–28.3 22.6 (10.4) 0.6–56.1 17.7 (8.7) 7.7–113.9 1214.2 (775.2) 10.2–3448.4 0.9 (0.2) 0.4–1 2.9 (0.8) 1.4–4.6 3.9 (0.9) 1.7–5 Bioregion # plots Hdom (m) BA (m2 ha−1) dg (cm) N (stems ha−1) SG ST DT Maritime 1044 17 (5.3) 3.5–33 28.1 (10.6) 0.7–78.2 26.1 (11.5) 8–111.5 755.6 (577) 5.1–3738.4 0.7 (0.1) 0.4–1 3.4 (1.1) 1.2–4.6 3 (0.8) 1.2–5 Mediterranean 374 12.4 (6) 3–27.6 22.8 (11.5) 0.4–52.6 19 (7) 7.6–81.3 998.8 (677.9) 10.2–2910.8 0.7 (0.2) 0.4–1 2.4 (0.8) 1.4–4.6 3.6 (1.1) 1.2–5 Mountainous 1121 17.8 (6) 2.5–36.5 31.5 (13.8) 0.4–82.7 25.6 (10) 7.6–95.5 808.2 (561.8) 5.1–3243.2 0.7 (0.1) 0.4–1 3.2 (1.4) 1.2–4.6 3.3 (0.9) 1.7–5 Transition 264 10 (4.1) 3–28.3 22.6 (10.4) 0.6–56.1 17.7 (8.7) 7.7–113.9 1214.2 (775.2) 10.2–3448.4 0.9 (0.2) 0.4–1 2.9 (0.8) 1.4–4.6 3.9 (0.9) 1.7–5 Average, standard deviation in parenthesis and range, minimum-maximum for Hdom: Dominant height, BA: basal area, dg: quadratic mean diameter, N: number of stems per hectare, SG: specific gravity according to Gutierrez Oliva and Plaza Pulgar (1967), ST: averaged shade tolerance, and DT: drought tolerance according to Niinimets and Valladares (2008). The dominant species and its respective SG, ST and DT values by bioregions are European beech (0.75, 4.56, 2.4) and Scots pine (0.50, 1.67, 4.34) in Maritime and Mountainous regions, Black pine (0.58, 2.1, 4.38) and Aleppo pine (0.55, 1.35, 4.97) in the Mediterranean region and Holm oak (1.00, 3.02, 4.72) and Faginea oak (0.95, 2.73, 3.02) in the Transition bioregion. Table 1 Descriptive stand variables, specific gravity and stress tolerances by bioregions in the study area. Bioregion # plots Hdom (m) BA (m2 ha−1) dg (cm) N (stems ha−1) SG ST DT Maritime 1044 17 (5.3) 3.5–33 28.1 (10.6) 0.7–78.2 26.1 (11.5) 8–111.5 755.6 (577) 5.1–3738.4 0.7 (0.1) 0.4–1 3.4 (1.1) 1.2–4.6 3 (0.8) 1.2–5 Mediterranean 374 12.4 (6) 3–27.6 22.8 (11.5) 0.4–52.6 19 (7) 7.6–81.3 998.8 (677.9) 10.2–2910.8 0.7 (0.2) 0.4–1 2.4 (0.8) 1.4–4.6 3.6 (1.1) 1.2–5 Mountainous 1121 17.8 (6) 2.5–36.5 31.5 (13.8) 0.4–82.7 25.6 (10) 7.6–95.5 808.2 (561.8) 5.1–3243.2 0.7 (0.1) 0.4–1 3.2 (1.4) 1.2–4.6 3.3 (0.9) 1.7–5 Transition 264 10 (4.1) 3–28.3 22.6 (10.4) 0.6–56.1 17.7 (8.7) 7.7–113.9 1214.2 (775.2) 10.2–3448.4 0.9 (0.2) 0.4–1 2.9 (0.8) 1.4–4.6 3.9 (0.9) 1.7–5 Bioregion # plots Hdom (m) BA (m2 ha−1) dg (cm) N (stems ha−1) SG ST DT Maritime 1044 17 (5.3) 3.5–33 28.1 (10.6) 0.7–78.2 26.1 (11.5) 8–111.5 755.6 (577) 5.1–3738.4 0.7 (0.1) 0.4–1 3.4 (1.1) 1.2–4.6 3 (0.8) 1.2–5 Mediterranean 374 12.4 (6) 3–27.6 22.8 (11.5) 0.4–52.6 19 (7) 7.6–81.3 998.8 (677.9) 10.2–2910.8 0.7 (0.2) 0.4–1 2.4 (0.8) 1.4–4.6 3.6 (1.1) 1.2–5 Mountainous 1121 17.8 (6) 2.5–36.5 31.5 (13.8) 0.4–82.7 25.6 (10) 7.6–95.5 808.2 (561.8) 5.1–3243.2 0.7 (0.1) 0.4–1 3.2 (1.4) 1.2–4.6 3.3 (0.9) 1.7–5 Transition 264 10 (4.1) 3–28.3 22.6 (10.4) 0.6–56.1 17.7 (8.7) 7.7–113.9 1214.2 (775.2) 10.2–3448.4 0.9 (0.2) 0.4–1 2.9 (0.8) 1.4–4.6 3.9 (0.9) 1.7–5 Average, standard deviation in parenthesis and range, minimum-maximum for Hdom: Dominant height, BA: basal area, dg: quadratic mean diameter, N: number of stems per hectare, SG: specific gravity according to Gutierrez Oliva and Plaza Pulgar (1967), ST: averaged shade tolerance, and DT: drought tolerance according to Niinimets and Valladares (2008). The dominant species and its respective SG, ST and DT values by bioregions are European beech (0.75, 4.56, 2.4) and Scots pine (0.50, 1.67, 4.34) in Maritime and Mountainous regions, Black pine (0.58, 2.1, 4.38) and Aleppo pine (0.55, 1.35, 4.97) in the Mediterranean region and Holm oak (1.00, 3.02, 4.72) and Faginea oak (0.95, 2.73, 3.02) in the Transition bioregion. Model structure The hypothesis was tested with a modified version of the model structure developed by Ducey et al. (2017). This model is based on the tree area ratio proposed by Chisman and Schumacher (1940) who defined the area occupied by the crown projection, i.e. the growing space requirement of a tree, as a quadratic function of stem diameter (equation (1)): Ai=c0+c1DBHi+c2DBHi2, (1) where Ai is the available growing space and DBHi is diameter at breast height of tree i, c0-c2 are parameters to be estimated. Summing all individual growing spaces results in the Tree Area Ratio (TAR, equation (2)) which represents the growing space occupied by trees in a stand. TAR=c0N+c1∑iDBHi+c2∑iDBHi2 (2) Although Chisman and Schumacher (1940) identified TAR with the potential crown space occupancy of a stand, other researchers have extended the interpretation to reflect a more abstract resource-based growing space requirement (sensuOliver and Larson, 1996). Setting TAR to 1 and estimating parameters c0,c1 and c2 with data from fully stocked stands, equation (2) can be used to predict the relative density of new stands. Ducey and Knapp (2010) took this approach but changed the growing space requirement to a function of the additive stand density index (Long and Daniel, 1990) combined with wood specific gravity following Woodall et al. (2005). More recently Ducey et al. (2017) modified the approach to accommodate covariates such as climate and functional traits like shade tolerance, e.g., Ai=(b0+∑j=1TbjTij)(DBHi25)1.6, (3) where Ai is the growing space of tree i, Tij is any plant trait of tree i (such as specific gravity or shade tolerance), DBHi is diameter at breast height of tree i and b0 through bj are parameters to be estimated. Summing over individual trees provides a relative density similar to TAR (equation (4)): RD=b0∑i=1n(DBHi25)1.6+∑j=1TbjTij∑i=1n(DBHi25)1.6. (4) Ducey et al. (2017) further extended the model to include climate variables as well as species functional traits. An appropriate fitting procedure is needed to estimate the parameters in equation (4). If, as in Chisman and Schumacher (1940), RD were set to equal 1 for all stands, and then least-squares regression were used to estimate the parameters, there would be an implicit assumption that the stands included in the data had an average stocking level corresponding to full stocking. Thus, it would be necessary to pre-screen the data to select only fully stocked stands. However, other fitting approaches do not need such pre-screening, as illustrated by Ducey and Knapp (2010) and Ducey et al. (2017) (see below in the Fitting procedure section). Stress tolerance and stress intensity Testing our hypothesis needs comparable coefficients and a reasonable way to modulate the stress tolerance as a function of the stress intensity. Our modification of Ducey et al.’s (2017) model includes a normalization of Niinemets and Valladares (2006) shade and drought tolerance rankings to make them comparable to the specific gravity values found in the dataset. We scaled DT and ST ranked values to a maximum of 1 by dividing the score value by the maximum value found in the dataset, renaming them as nDT and nST, respectively. In so doing, we can assess the relative contribution of each tolerance to the maximum carrying capacity of the stand as all three traits have a maximum value of 1 (e.g. maximum SG is 1 g cm−3 for holm oak, Quercus ilex L.). Instead of using continuous climatic variables to assess the climate effect on stress tolerance we built multipliers which reduced the normalized stress tolerance in harsh conditions, hereafter referred as stress tolerance modulators, allowing us to order species tolerance from 0 (no stress and the tolerance value does not influence the maximum density) to the maximum tolerance, i.e. the maximum normalized tolerance value nDT or nST. The contribution of DT to predict the maximum density is regulated by a multiplier that maximizes the tolerance from lower to higher drought stress levels: dI=kmax(k), (5) where dI is the drought intensity stress modulator of nDT, k is the drought intensity, measured for each plot as the area confined between the temperature and precipitation curves in a Walter-Lieth’s climodiagram, max(k) is the maximum drought intensity for the study area, this modulator ranges from 0, i.e. sites without drought period where the contribution of drought tolerance is negligible (k = 0), to 1, i.e. sites with maximum drought intensity (k = max(k)). Analogously, shade tolerance is modulated by minimum temperature, which is correlated with the length of the growing season (Laanisto and Niinemets, 2015; Ducey et al., 2017): cI=1−Tmin−minTminmaxTmin−minTmin, (6) where cI is the cold-stress modulator of nST, Tmin is the minimum temperature of the coldest month, minTmin and maxTmin are the minimum and the maximum of the minimum temperature of the coldest month, respectively, this modulator ranges from 1 (minimum minT, associated to cold sites where shade tolerance is maximum for the species) to 0 (warmer sites). Although specific gravity might vary, among others, with altitude and dry conditions (Chave et al., 2006; Pompa-García and Venegas-González, 2016) or exposure to stressors, like flooding (Lawson et al., 2015) we assume, according to our working hypothesis 1, that bending stress is always occurring across the species distribution as it determines the maximum amount of foliage that a stand can support (Dean and Baldwin, 1996). In addition, there is no significant relationship between specific gravity and site aridity measured by de Martonne Index in the study area (see Supplementary Figure 1). The trade-offs and synergies of stress tolerance modulators were analysed with dispersion plots of two contrasting broadleaves (Fagus sylvatica L. and Quercus ilex L.) and conifers (Pinus sylvestris L. and P. halepensis Mill.). A species with great variability in both multipliers, i.e. large dispersion in both axes, would indicate a large suitable habitat where there could be a shift in the relative importance of shade or drought tolerance and even polytolerance, whereas low variability in either of both multipliers would indicate absence of polytolerance and that stress tolerance to a single stressor is more important for the species to occupy the site. Hypothesis testing We started the fitting procedure from the base model with SG as a single variable as in Ducey and Knapp (2010), equation (4) (hypothesis 1), and adding one more normalized stress tolerance alternatively (nST and nDT) multiplied by its corresponding stress modulator (cI and dI, respectively). The subsequent models used for the rest of hypotheses are as follows: Testing hypothesis 2 expand equation (4) with a term including the tolerance to shade RD=a0∑iEFi(DBHi25)1.6+a1∑iEFiSGi(DBHi25)1.6+a2∑iEFi⋅(cI⋅nSTi)⋅(DBHi25)1.6, (7) where RD is relative density, EFi is the expansion factor applicable to each individual according to its position in the concentric plot design of the NFI, DBHi is the diameter of individual, SGi is the species-specific wood gravity for each individual, nSTi is the species-specific normalized shade tolerance for each individual, cI is the cold temperature stress modulator, and a0, a1 and a2 are parameters to be estimated. Hypothesis 3 was tested by adding a term for drought tolerance to equation (4): RD=a0∑iEFi(DBHi25)1.6+a1∑iEFiSGi(DBHi25)1.6+a3∑iEFi⋅(dI⋅nDTi)(DBHi25)1.6, (8) where Di is the drought stress modulator, nDTi is the species-specific normalized drought tolerance and the rest of symbols are as in equation (7). Alternatively, we tested the simultaneous relative importance of tolerance to all stressors: RD=a0∑iEFi(DBHi25)1.6+a1∑iEFiSGi(DBHi25)1.6+a2∑iEFi⋅(cI⋅nSTi)⋅(DBHi25)1.6+a3∑iEFi⋅(dI⋅nDTi)(DBHi25)1.6, (9) where symbols were previously defined. Model specification described in equations (5) and (6) does not allow for simultaneous drought and shade tolerance values in neither cold nor warm sites (i.e. polytolerance). We also fitted the complementary value in the cold index (1-cI) to allow for testing shade–drought tolerance in warmer sites (equation (10)): RD=a0∑iEFi(DBHi25)1.6+a1∑iEFiSGi(DBHi25)1.6+a2∑iEFi⋅((1−cI)⋅nSTi)⋅(DBHi25)1.6+a3∑iEFi⋅(dI⋅nDTi)(DBHi25)1.6, (10) where all symbols were previously defined. Fitting procedure Traditional fitting procedures of maximum stocking need selection of plots at full stocking levels. However, screening data from national forest inventory plots is not efficient and rather subjective as the stocking level for RD = 1 (full site occupancy) is not truly known and there is evidence that plots in the data set do not correspond to RD = 1 (Ducey and Knapp, 2010), especially in large areas where environmental conditions and species traits can vary affecting the degree of occupancy. Plantations can be an exception and consequently they were removed from the data set, representing less than 1 per cent of the original dataset. Fitting models derived from equation (4) by ordinary least squares would require an independent and normal error term, and the regression procedure would select coefficients such that the mean value of the plots included in the dataset corresponded to RD = 1. However, our purpose here is to find a combination of coefficients such that the calculated RD value would fall below RD = 1 for nearly all plots; in other words, we seek a maximum density relationship at RD = 1 that would fall above all observations if the observations were free of measurement errors and other anomalies. To avoid subjective pre-screening procedures when using National Forest Inventory datasets, a default value of RD = 1 to all plots is hypothesized, and the models (equations (7) through (10)) are fit by quantile regression instead of ordinary least squares. The hypothetical assignment of RD = 1 to all plots does not mean that the actual density of all plots is at the maximum. Instead, it is giving a reference value of RD that should be consistent across all plots (no matter what their tree list, species traits, and climate might be). In such case, quantile regression allows the actual density of individual plots to vary freely, but it does enforce a model fit so that the reference RD = 1 actually falls at an extreme quantile of the distribution of RD values for individual plots, rather than the mean. By incorporating both species traits and climate variables and predictors, quantile regression ensures that the reference RD = 1, taken as a conditional quantile, remains consistently near the upper end of the distribution. Quantile regression fits a linear model through any value in data as long as the modeller selects a quantile. A quantile is a proportion, q∈[0,1] splitting the data into two proportions where the q observations are below the regression line and the 1-q are above. Maximum density needs a high quantile in order to assure that plots are at maximum density when on the modelled line; occasional observations may fall above the line due to measurement errors or unusual spatial configurations of trees relative to the plot boundaries. We tested four q values leaving 0.5 per cent, 1 per cent, 2.5 per cent and 5 per cent of observations above the maximum density boundary, i.e. we fitted quantiles 0.995, 0.99, 0.975 and 0.95. By setting RD = 1 to all plots and using the quantile regression we are removing all subjectivity associated to plot selection. Quantile regression has been successfully applied to model maximum boundary lines (Zhang et al., 2005; Ducey and Knapp, 2010; Vospernik and Sterba, 2014; Ducey et al., 2017). Schwarz’s Bayesian Information Criteria (BIC) within each quantile was used to assess fitting performance: BIC=−2⋅ln(L)−k⋅ln(n), (11) where L is the log-likelihood of the data, k is the number of parameters and n the number of observations. The lower BIC the better the fit. The significance level for parameter estimation was set to 0.05. All models were fitted using quantreg package of R (Koenker, 2016; R Core Team, 2016). Evaluation of maximum SDI for pure and mixed stands Maximum density predictions were evaluated graphically across environmental stress gradients (cold temperature and drought) between and within functional groups and species composition (broadleaves vs conifers and pure vs mixed-species stands). The maximum stand density relationship for two important commercial species, Scots pine (Pinus sylvestris L.) and European beech (Fagus sylvatica L.), has been previously assessed in the study area by Condés et al. (2017, 2013). Their maximum SDI was used as benchmark against to discuss our model’s output. In order to contrast results against climatic regions we also included comparisons of two Mediterranean species Aleppo pine (Pinus halepensis Mill.) and Holm oak (Quercus ilex L.) Additionally, we compared the SDI predicted by our model with information of SDI retrieved in the literature for species also found in our database: Pedunculate and Sessile oak (Quercus robur and Quercus petraea), black pine (Pinus nigra Arn.) and Silver fir (Abies alba L.) using values published by Pretzsch and Biber (2005) and Vospernik and Sterba (2014). Results Relationship and variability between normalized stress tolerance values Simultaneous variation in tolerance to drought and cold would indicate that a species or group of species can be found in a wider environmental range. Differences in polytolerance to shade and drought were detected by functional groups (Figure 2a and b). Plots of broadleaves growing in pure and mixed stands showed the highest variability for simultaneous normalized drought and shade tolerance whereas conifers did not show polytolerance either in pure or in mixed stands. Figure 2 View largeDownload slide Relationship between drought and shade tolerance ranking modulated by drought and cold temperature stress intensity (dI and cI). Large variability in both axes would indicate polytolerance. Figure 2 View largeDownload slide Relationship between drought and shade tolerance ranking modulated by drought and cold temperature stress intensity (dI and cI). Large variability in both axes would indicate polytolerance. A graphical analysis of stress modulators for two contrasting species per functional group showed that Scots pine and European beech are located mostly on drought free sites but with a large variability of minimum temperatures indicating that both species are not drought tolerant, whereas Holm oak (Quercus ilex L.) and Aleppo pine (Pinus halepensis Mill.) can be found over both a large range of drought intensity values and minimum temperature range (Figure 3a and b). Figure 3 View largeDownload slide (a) Relationship between cold and drought multipliers for broadleaves in Navarra (Spain). Large point dispersion in both axes would indicate simultaneous tolerance for cold and drought stress. Filled circles for Fagus sylvatica, void circles for Quercus ilex. (b) Relationship between cold and drought multipliers for conifers in Navarra (Spain). Large point dispersion in both axes would indicate simultaneous tolerance for cold and drought stress. Filled circles for Pinus sylvestris, void circles for Pinus halepensis. Figure 3 View largeDownload slide (a) Relationship between cold and drought multipliers for broadleaves in Navarra (Spain). Large point dispersion in both axes would indicate simultaneous tolerance for cold and drought stress. Filled circles for Fagus sylvatica, void circles for Quercus ilex. (b) Relationship between cold and drought multipliers for conifers in Navarra (Spain). Large point dispersion in both axes would indicate simultaneous tolerance for cold and drought stress. Filled circles for Pinus sylvestris, void circles for Pinus halepensis. Best model selection SG as the only explanatory variable (Table 2) was significant at all quantiles (equation (4), hypothesis 1 accepted) leading to an implied maximum density for Scots pine ranging from 1427 stems ha−1 at quantile 0.995 to 1121 stems ha−1 at quantile 0.95. European beech showed a maximum density between 1233 stems ha−1 and 945 stems ha−1 for the same quantiles, respectively. The addition of shade tolerance alone led to non-significant parameters at 0.05 significant level (equation (7), hypothesis 2 rejected), whereas drought tolerance alone led to significant models at all quantiles (equation (8), hypothesis 3 accepted). The full model with ST acting in cold sites and DT in warm sites (equation (9)) led to non-significant parameters for ST at all quantiles. The best fitted model in terms of BIC includes both normalized stress tolerance values modulated by temperature and drought intensity in warmer sites (equation (10)). This model at quantile 0.975 will be used for further SDI predictions as it represents the best fitted upper boundary line. Table 2 Estimates of relative density models including functional traits (SG, specific gravity; DT, drought tolerance; ST, shade tolerance) and abiotic stress modulators (iC cold temperature modulator and iD drought modulator). Model Variable q0.995 q0.99 q0.975 q0.95 Estimate (s.e.) BIC Estimate (s.e.) BIC Estimate (s.e.) BIC Estimate (s.e.) BIC (4) Intercept 0.00048 (0.00009) 3364.3 0.00043 (0.0001) 3305.9 0.00044 (0.00011) 3206.57 0.00051 (0.0001) 3113.18 SG 0.00044 (0.00013) 0.00059 (0.00016) 0.00068 (0.00015) 0.00073 (0.00013) (7) Intercept 0.00048 (0.00009) 0.0004 (0.0001) 0.00044 (0.00012) 0.00051 (0.0001) SG 0.00044 (0.00013) 0.00061 (0.00017) 0.00068 (0.00015) 0.00078 (0.00014) iC x ST n.s. n.s. n.s. n.s. (8) Intercept 0.00048 (0.00008) 3317.44 0.00043 (0.0001) 3273.17 0.00044 (0.00011) 3171.97 0.00051 (0.00009) 3077.7 SG 0.00044 (0.00013) 0.00059 (0.00016) 0.00068 (0.00016) 0.00071 (0.00013) iD x DT 0.00129 (0.00033) 0.00126 (0.00022) 0.0012 (0.00022) 0.00121 (0.00028) (9) Intercept 0.00048 (0.00008) 0.0004 (0.0001) 0.00014 (0.00012) 0.00051 (0.00009) SG 0.0044 (0.00011) 0.00061 (0.00018) 0.00068 (0.00015) 0.00076 (0.00013) iC x ST n.s. n.s. n.s. n.s. iD x DT 0.00129 (0.00039) 0.00127 (0.00023) 0.0012 (0.00025) 0.00115 (0.00026) (10) Intercept 0.00049 (0.00008) 0.00051 (0.00011) 0.00052 (0.0001) 3116.91 0.00054 (0.00009) 3035.1 SG n.s. n.s. 0.00033 (0.00016) 0.00045 (0.00014) (1−iC) x ST 0.00029 (0.00013) n.s. 0.00046 (0.00015) 0.00041 (0.00011) iD x DT 0.0013 (0.00021) 0.00128 (0.00023) 0.00122 (0.0002) 0.00115 (0.00029) Model Variable q0.995 q0.99 q0.975 q0.95 Estimate (s.e.) BIC Estimate (s.e.) BIC Estimate (s.e.) BIC Estimate (s.e.) BIC (4) Intercept 0.00048 (0.00009) 3364.3 0.00043 (0.0001) 3305.9 0.00044 (0.00011) 3206.57 0.00051 (0.0001) 3113.18 SG 0.00044 (0.00013) 0.00059 (0.00016) 0.00068 (0.00015) 0.00073 (0.00013) (7) Intercept 0.00048 (0.00009) 0.0004 (0.0001) 0.00044 (0.00012) 0.00051 (0.0001) SG 0.00044 (0.00013) 0.00061 (0.00017) 0.00068 (0.00015) 0.00078 (0.00014) iC x ST n.s. n.s. n.s. n.s. (8) Intercept 0.00048 (0.00008) 3317.44 0.00043 (0.0001) 3273.17 0.00044 (0.00011) 3171.97 0.00051 (0.00009) 3077.7 SG 0.00044 (0.00013) 0.00059 (0.00016) 0.00068 (0.00016) 0.00071 (0.00013) iD x DT 0.00129 (0.00033) 0.00126 (0.00022) 0.0012 (0.00022) 0.00121 (0.00028) (9) Intercept 0.00048 (0.00008) 0.0004 (0.0001) 0.00014 (0.00012) 0.00051 (0.00009) SG 0.0044 (0.00011) 0.00061 (0.00018) 0.00068 (0.00015) 0.00076 (0.00013) iC x ST n.s. n.s. n.s. n.s. iD x DT 0.00129 (0.00039) 0.00127 (0.00023) 0.0012 (0.00025) 0.00115 (0.00026) (10) Intercept 0.00049 (0.00008) 0.00051 (0.00011) 0.00052 (0.0001) 3116.91 0.00054 (0.00009) 3035.1 SG n.s. n.s. 0.00033 (0.00016) 0.00045 (0.00014) (1−iC) x ST 0.00029 (0.00013) n.s. 0.00046 (0.00015) 0.00041 (0.00011) iD x DT 0.0013 (0.00021) 0.00128 (0.00023) 0.00122 (0.0002) 0.00115 (0.00029) Table 2 Estimates of relative density models including functional traits (SG, specific gravity; DT, drought tolerance; ST, shade tolerance) and abiotic stress modulators (iC cold temperature modulator and iD drought modulator). Model Variable q0.995 q0.99 q0.975 q0.95 Estimate (s.e.) BIC Estimate (s.e.) BIC Estimate (s.e.) BIC Estimate (s.e.) BIC (4) Intercept 0.00048 (0.00009) 3364.3 0.00043 (0.0001) 3305.9 0.00044 (0.00011) 3206.57 0.00051 (0.0001) 3113.18 SG 0.00044 (0.00013) 0.00059 (0.00016) 0.00068 (0.00015) 0.00073 (0.00013) (7) Intercept 0.00048 (0.00009) 0.0004 (0.0001) 0.00044 (0.00012) 0.00051 (0.0001) SG 0.00044 (0.00013) 0.00061 (0.00017) 0.00068 (0.00015) 0.00078 (0.00014) iC x ST n.s. n.s. n.s. n.s. (8) Intercept 0.00048 (0.00008) 3317.44 0.00043 (0.0001) 3273.17 0.00044 (0.00011) 3171.97 0.00051 (0.00009) 3077.7 SG 0.00044 (0.00013) 0.00059 (0.00016) 0.00068 (0.00016) 0.00071 (0.00013) iD x DT 0.00129 (0.00033) 0.00126 (0.00022) 0.0012 (0.00022) 0.00121 (0.00028) (9) Intercept 0.00048 (0.00008) 0.0004 (0.0001) 0.00014 (0.00012) 0.00051 (0.00009) SG 0.0044 (0.00011) 0.00061 (0.00018) 0.00068 (0.00015) 0.00076 (0.00013) iC x ST n.s. n.s. n.s. n.s. iD x DT 0.00129 (0.00039) 0.00127 (0.00023) 0.0012 (0.00025) 0.00115 (0.00026) (10) Intercept 0.00049 (0.00008) 0.00051 (0.00011) 0.00052 (0.0001) 3116.91 0.00054 (0.00009) 3035.1 SG n.s. n.s. 0.00033 (0.00016) 0.00045 (0.00014) (1−iC) x ST 0.00029 (0.00013) n.s. 0.00046 (0.00015) 0.00041 (0.00011) iD x DT 0.0013 (0.00021) 0.00128 (0.00023) 0.00122 (0.0002) 0.00115 (0.00029) Model Variable q0.995 q0.99 q0.975 q0.95 Estimate (s.e.) BIC Estimate (s.e.) BIC Estimate (s.e.) BIC Estimate (s.e.) BIC (4) Intercept 0.00048 (0.00009) 3364.3 0.00043 (0.0001) 3305.9 0.00044 (0.00011) 3206.57 0.00051 (0.0001) 3113.18 SG 0.00044 (0.00013) 0.00059 (0.00016) 0.00068 (0.00015) 0.00073 (0.00013) (7) Intercept 0.00048 (0.00009) 0.0004 (0.0001) 0.00044 (0.00012) 0.00051 (0.0001) SG 0.00044 (0.00013) 0.00061 (0.00017) 0.00068 (0.00015) 0.00078 (0.00014) iC x ST n.s. n.s. n.s. n.s. (8) Intercept 0.00048 (0.00008) 3317.44 0.00043 (0.0001) 3273.17 0.00044 (0.00011) 3171.97 0.00051 (0.00009) 3077.7 SG 0.00044 (0.00013) 0.00059 (0.00016) 0.00068 (0.00016) 0.00071 (0.00013) iD x DT 0.00129 (0.00033) 0.00126 (0.00022) 0.0012 (0.00022) 0.00121 (0.00028) (9) Intercept 0.00048 (0.00008) 0.0004 (0.0001) 0.00014 (0.00012) 0.00051 (0.00009) SG 0.0044 (0.00011) 0.00061 (0.00018) 0.00068 (0.00015) 0.00076 (0.00013) iC x ST n.s. n.s. n.s. n.s. iD x DT 0.00129 (0.00039) 0.00127 (0.00023) 0.0012 (0.00025) 0.00115 (0.00026) (10) Intercept 0.00049 (0.00008) 0.00051 (0.00011) 0.00052 (0.0001) 3116.91 0.00054 (0.00009) 3035.1 SG n.s. n.s. 0.00033 (0.00016) 0.00045 (0.00014) (1−iC) x ST 0.00029 (0.00013) n.s. 0.00046 (0.00015) 0.00041 (0.00011) iD x DT 0.0013 (0.00021) 0.00128 (0.00023) 0.00122 (0.0002) 0.00115 (0.00029) According to the best model, the relative contribution of shade and drought tolerance to explain the maximum density for a given species increases in warm sites. Specific wood gravity and its relation to bending stress contributes more to explain maximum density in colder sites. Predicted maximum SDI values Differences were found in maximum density between phylogenetic groups. Gymnosperms showed higher maximum number of stems per unit area at a given diameter than angiosperms. Within each phylogenetic group drought tolerant species had lower density than drought intolerant species whereas shade intolerant species had higher density. The maximum density is reduced across the environmental stress gradient with the abiotic stress factor changing from the south where drought is more intense to north, where cold temperatures is the limiting factor. The maximum SDI for Scots pine ranges from 1244 to 1416 stems ha−1, whereas European beech maximum density ranged from 826 to 1191 stems ha−1 (Figure 4a). Drought tolerant species needed fewer stems to full occupy the stand. Aleppo pine ranged from 473 to 1225 stems ha−1 and Holm oak ranged from 725 to 946 stems ha−1 (Figure 4b). Figure 4 View largeDownload slide (a) Predicted maximum stand density index (stems ha−1) across temperature stress gradient in Navarra, Spain, for cold-tolerant/drought intolerant species. Grey circles for Pinus sylvestris and black circles for Fagus sylvatica. (b) Predicted maximum stand density index (stems ha−1) across drought stress gradient in Navarra, Spain, for cold-intolerant/drought tolerant species. Grey circles for Pinus halepensis and black circles for Quercus ilex. Figure 4 View largeDownload slide (a) Predicted maximum stand density index (stems ha−1) across temperature stress gradient in Navarra, Spain, for cold-tolerant/drought intolerant species. Grey circles for Pinus sylvestris and black circles for Fagus sylvatica. (b) Predicted maximum stand density index (stems ha−1) across drought stress gradient in Navarra, Spain, for cold-intolerant/drought tolerant species. Grey circles for Pinus halepensis and black circles for Quercus ilex. Discussion Maximum density is important to our understanding of species stand dynamics as it indicates the maximum occupancy or full-resource use by trees. This relationship has been described as the link between forestry and ecology (Jack and Long, 1996; Long et al., 2004) and it relies in the constancy of leaf area index (LAI) during self-thinning (Long and Smith, 1984). The assumption of constancy of LAI has helped to build the hypothesis that tolerance to bending stress is driven by structural traits like specific wood density. Woodall et al. (2005) demonstrated that under this assumption trees with lower wood density would need fewer individuals to support the same amount of foliage than species with higher wood density. This hypothesis could be behind the evidence that Scots pine and beech mixed stands in Europe have greater densities than pure stands (Pretzsch and Schütze, 2015; Pretzsch et al., 2015) as this type of mixture also has lower wood density (Zeller et al., 2017). So far, this hypothesis has been tested for North America species and we present here the first test of this hypothesis for European tree species expanding the validity of the stress-bending stress control over stand density. However, whether LAI remains constant during self-thinning has been put into question as foliar mass development is age-dependent indicating that the leaf mass that a stand can support is a function of stand development (Holdaway et al., 2008). These authors pointed that the change in leaf mass can be attributable to changes in crown vertical structure. On the other hand, Pretzsch and Mette (2008) found constant LAI in self-thinned stands of Norway spruce and beech but variable specific leaf area (SLA) with increasing diameter attributed to losses of shade leaves from suppressed dead trees. Self-thinning implies redistribution of the same amount of foliage in less and larger individuals (Long et al. 2004). However, the non-constancy of LAI is more evident in old stands where large aged trees cannot occupied the space available after mortality or thinning events leading to departures from a linear size-density trajectory (Long and Vacchiano, 2014). In general, the older the individual tree the less the foliage in high dense stands (Long et al. 2004). To reconcile the LAI constancy hypothesis with observations of age-dependent LAI, self-thinning boundary line modelling must combine the abiotic stress intensity and tolerance approach, as done in our model, and the existing relationship between crown plasticity and shade tolerance (Pretzsch, 2014; Jucker et al., 2015) or the competitive advantage of species with higher crown plasticity (Vincent and Harja, 2008). Equation (10) predicts a change in the relative contribution of stress tolerance from specific wood gravity in cold sites to increasing contribution of drought and shade tolerance in dry sites. This is in line with the finding of Condés et al. (2017) who suggested that plasticity in hydraulic structure outweighs mechanical control over stand density in arid environments. However, tree hydraulics might not be only relevant in dry areas as water transport has been also proposed to explain the inverse relationship between wood density and leaf mass in Neotropical forests (Wright et al., 2007). Abiotic stress tolerance is an important driver of woody plant distribution and productivity, as negative changes of environmental conditions may reduce the realized niche and potential biomass growth of a species (Laanisto and Niinemets, 2015). Polytolerance to stress is considered rare in nature leading to a trade-off between shade and drought tolerance (Niinemets and Valladares, 2006). Equation (9) in Table 2 would have yielded a lack of simultaneous effect of shade and drought tolerance in maximum density at both ends of a cold-warm spectrum (Supplementary Figure 2a). However, the lack of statistical significance for shade tolerance in equation (9) and the significance of its complementary value in equation (10) indicates that the relative importance of shade and drought tolerance increases in warmer sites (Supplementary Figure 2b). Tolerance to both stressors result in a reduction of the number of individuals needed to make full use of resources. However, species growing in warm sites are usually more drought tolerant than shade tolerant and the predicted reduction is mainly driven by drought stress, as it was reflected by the relative contribution of drought tolerance which is 2.7 higher than shade tolerance (Table 2). In contrast, species growing in cold sites showed an increasing effect of tolerance to bending stress on the number of individuals as compared to tolerance to shade and drought. Cold sites are expected to have a shorter growing season which has been pointed out to be a main reason for shade–drought trade-offs (Valladares and Niinemets, 2008; Laanisto and Niinemets, 2015). Equation (10) does not predict the presence or absence of polytolerance but it can be used to make inferences on maximum density about the effect of observed shade–drought trade-offs or polytolerance in a cold-warm spectrum. Warm sites would support less trees per unit area irrespective of tolerance ranking whereas a suite of polytolerant species would need fewer individuals than species with stronger trade-offs (Supplementary Figure 3). Linking this finding with the fact that polytolerance is associated with a higher probability of survival in harsh conditions, as suggested by Laanisto and Niinemets (2015), would have to be assessed in terms of constancy or stress tolerance dependence of the slope parameter in the density-size relationship. The basic model used in our approach keeps a constant slope in the density-size relationship (equation (4)). The constancy of slope in forestry and ecology has been a subject of intense debate (Lonsdale, 1990; Zeide, 2004, 1987) and findings suggest that the 1.6 constant slope is much more an exception than a rule (del Rı́o et al., 2001; Pretzsch and Biber, 2005). The mathematical computation of a variable slope according to environmental filters or stress tolerances in our model is challenging but it would improve our understanding of the mechanisms behind site occupancy. Nevertheless, our SDI predictions using stress intensity and tolerance are realistic and do not dramatically depart from predictions using a variable slope (del Rı́o et al., 2001; Pretzsch and Biber, 2005; Vospernik and Sterba, 2014; Condés et al., 2017). SDI for Scots pine in our study, 1307 stems ha−1, is consistent with 1297 stems ha−1 found in Condés et al. (2013) for the same region and slightly less than the 1442 stems ha−1 found by Río et al. (2008) for Spain. For European beech our model predicted less individuals, 966 stems ha−1, as compared with 1042 stems ha−1 found by Condés et al. (2013), nevertheless the range of our predictions embraced that of Condés’ as it ranged from 826 in colder sites to 1191 stems ha−1 in warmer sites. Our model predictions are closer to other estimations for beech in Europe (Figure 5), whereas SDI values for Silver fir (1134 stems ha−1) and Black pine (1246 stems ha−1) in Navarra are sightless less from those found in Austria, 1278 and 1424 stems ha−1, respectively (Vospernik and Sterba, 2014). Figure 5 View largeDownload slide Mean SDI (bars), standard error (whiskers), maximum SDI (triangles) found in the literature search for the species in the x-axis. Black points are predictions in Navarra using model 4 in Table 2. Figure 5 View largeDownload slide Mean SDI (bars), standard error (whiskers), maximum SDI (triangles) found in the literature search for the species in the x-axis. Black points are predictions in Navarra using model 4 in Table 2. Stress tolerance to shade is an important driver in shaping species distribution and community composition whereas structural traits, like wood density predict species’ competitive response (Kunstler et al., 2016). However, tolerance is not constant for a given species and it varies across genotypes and environmental conditions (Barnes et al., 1997). For example, shade tolerance is affected by nutrient and water availability (Valladares and Niinemets, 2008 and reference therein) and it decreases for Douglas-fir, western redcedar and hemlock with increasing soil moisture. Stress tolerance differs both with species ontogeny being more important in juvenile stages than in mature forests, and might differ during the same stage within a gradient of stress and species with a large distribution area. For example, Scots pine does not show the same shade tolerance across its distribution in Europe, being more shade intolerant in Central Europe than in the South. This variation in stress tolerance allows species to survive and grow in a wide range of light conditions. This gives a clue that scoring species according to stress tolerance is not an absolute value and the tolerance might be modulated by climatic factors as tested by Ducey et al. (2017) and normalized in this paper. Our model cannot be used to test the over- or under-density of mixed vs pure stands as the contribution of stress tolerance to maximum density is averaged across species, i.e. the resulting density cannot be higher that the individual contribution of each species. This is a caveat as long as several studies have shown higher density in mixed stands than in single-species ones (Pretzsch et al., 2015; Pretzsch and Biber, 2016) in association to changes in specific wood density in mixed stands (Zeller et al., 2017). A practical use of the model presented here is the substitution of stand density measures in empirical growth and yield model to include stress tolerance in the prediction of timber production. In so doing, the stand level management models will capture differences in yield production according to differences in stress tolerances. Silviculture aims to maintain density below the maximum density to keep stocking between the onset of competition and the onset of the density-dependent mortality where growth is constrained by resource competition and production loss due to mortality is not economically important. When management objectives promote individual tree growth the density is kept near the lower limit but if higher stand level yield is expected then the density is kept near the upper limit. These limits are set between 35 and 65 per cent of maximum density based on the seminal paper by Long (1985) and have been used in silvicultural prescriptions (Valbuena et al., 2008), however these values remain empirical and a mechanistic explanation is still lacking. Although our model does not provide mechanistic information about the exact limits it helps to accommodate the density taking into account the abiotic stress intensity and the species tolerance. For example, in the case of managed Scots pine stands in the study area the SDI limits should be set between 868 and 467 stems ha−1 in cold sites and between 544 and 293 stems ha−1 in the warmest sites. Conclusions We built a normalized approach to explain the maximum density needed to make full use of resources in terms of variability of the stress tolerance ranking according to Niinemets and Valladares (2006) modulated by the intensity of drought and low temperature stress. The model is based on the assumption that bending stress is concomitant to all sites in the study area assuming a constant LAI during self-thinning. A basic model of stocking degree inversely related to specific wood density is found significant in the study area showing an inverse relationship (hypothesis 1 accepted). The addition of normalized tolerance shade showed to be non-significant (hypothesis 2 rejected) whereas the combination of drought tolerance reduces the number of individuals that a site can support being this reduction stronger in warm sites (hypothesis 3 accepted). A species or group of species with high capacity to tolerate several abiotic stresses simultaneously, i.e. polytolerance, need far less individuals to fully occupy a site than species showing tolerance trade-offs. Linking polytolerance/trade-offs with a mechanistic explanation of survival at the stand level needs further investigation that should include variation due to stress tolerance in the slope parameter of the density-size relationship. Supplementary data Supplementary data are available at Forestry online. Acknowledgements The authors thank the anonymous reviewers and handling Editor for their constructive comments. Olivia L. Bartlett provided additional helpful suggestions on the manuscript. Funding The Spanish Minister of Economy and Competitiveness (grant number: AGL2014–51964-C2–2-R); the New Hampshire Agricultural Experiment Station (Scientific Contribution Number 2768); the USDA National Institute of Food and Agriculture McIntire-Stennis (Project Accession Number 1007007); and Short Term Scientific Mission provided by COST Action EuMIXFOR (COST Action number FP1206) to A.B.-O. Conflict of interest statement None declared. References Assmann , E . 1971 The Principles of Forest Yield Study . Pergamon press . Barnes , B.V. , Zak , D.R. , Denton , S.R. and Spurr , S.H. 1997 Forest Ecology . 4th ed . John Wiley and Sons . Burkhart , H.E. and Tomé , M. 2012 Modeling Forest Trees and Stands . Springer . Google Scholar CrossRef Search ADS Chave , J. , Muller-Landau , H.C. , Baker , T.R. , Easdlae , T.A. , Ter Steege , H. , et al. 2006 Regional and phylogenetic variation of wood denisty across 2456 neotropical tree species . Ecol. Appl. 16 , 2356 – 2367 . doi:10.3159/1095-5674(2007)134[301:SCROHD]2.0.CO;2 . Google Scholar CrossRef Search ADS PubMed Chisman , H.H. , and Schumacher , F.X. 1940 On the tree-area ratio and certain of its applications . J. For. 38 ( 4 ), 311 – 317 . Condés , S. , Del , M. and Sterba , H. 2013 Mixing effect on volume growth of Fagus sylvatica and Pinus sylvestris is modulated by stand density . For. Ecol. Manag. 292 , 86 – 95 . Google Scholar CrossRef Search ADS Condés , S. , Vallet , P. , Bielak , K. , Bravo-Oviedo , A. , Coll , L. , Ducey , M.J. , et al. 2017 Climate influences on the maximum size-density relationship in Scots pine (Pinus sylvestris L.) and European beech (Fagus sylvatica L.) stands . For. Ecol. Manag. 385 , 295 – 307 . doi:10.1016/j.foreco.2016.10.059 . Google Scholar CrossRef Search ADS Curtis , R.O. 1971 A tree area power function and related stand density measures for Douglas-fir . For. Sci. 17 , 146 – 159 . de Martonne , E. 1926 Une Nouvelle Fonction CIimatologique: L’Indice d’Aridite . La Météorol. 2 , 449 – 458 . Dean , T.J. and Baldwin , V.C. 1996 The relationship between Reineke’s stand-density index and physical stem mechanics . For. Ecol. Manag. 81 , 25 – 34 . doi:10.1016/0378-1127(95)03666-0 . Google Scholar CrossRef Search ADS del Rı́o , M. , Montero , G. and Bravo , F. 2001 Analysis of diameter–density relationships and self-thinning in non-thinned even-aged Scots pine stands . For. Ecol. Manag. 142 , 79 – 87 . doi:10.1016/S0378-1127(00)00341-8 . Google Scholar CrossRef Search ADS Ducey , M.J. and Knapp , R.A. 2010 A stand density index for complex mixed species forests in the northeastern United States . For. Ecol. Manag. 260 , 1613 – 1622 . doi:10.1016/j.foreco.2010.08.014 . Google Scholar CrossRef Search ADS Ducey , M.J. , Woodall , C.W. and Bravo-Oviedo , A. 2017 Climate and species functional traits influence maximum live tree stocking in the Lake States, USA . For. Ecol. Manag. 386 , 51 – 61 . doi:10.1016/j.foreco.2016.12.007 . Google Scholar CrossRef Search ADS Gonzalo-Jiménez , J. 2010 Diagnosis fitoclimática de la España Peninsular: hacia un modelo de clasificación funcional de la vegetación y de los ecosistemas peninsulares españoles. Edited by Organismo Autónomo de Parques Nacionales Gutiérrez , A. , and Plaza , F. 1967 Características fisico-mecánicas de las maderas Españolas. Madrid: Ministerio de Agricultura. Dirección General de Montes, Caza y Pesca Fluvial. IFIE . Holdaway , R.J. , Allen , R.B. , Clinton , P.W. , Davis , M.R. and Coomes , D.A. 2008 Intraspecific changes in forest canopy allometries during self-thinning . Funct. Ecol. 22 , 460 – 469 . doi:10.1111/j.1365-2435.2008.01388.x . Google Scholar CrossRef Search ADS Jack , S.B. and Long , J.N. 1996 Linkages between silviculture and ecology: an analysis of density management diagrams . For. Ecol. Manag. 86 , 205 – 220 . doi:10.1016/S0378-1127(96)03770-X . Google Scholar CrossRef Search ADS Jucker , T. , Bouriaud , O. and Coomes , D.A. 2015 Crown plasticity enables trees to optimize canopy packing in mixed-species forests . Funct. Ecol. 29 , 1078 – 1086 . doi:10.1111/1365-2435.12428 . Google Scholar CrossRef Search ADS Koenker , R. , 2016 . quantreg: Quantile Regression. R package version 5.29. Kunstler , G. , Falster , D. , Coomes , D.A. , Hui , F. , Kooyman , R.M. , Laughlin , D.C. , et al. 2016 Plant functional traits have globally consistent effects on competition . Nature 529 , 204 – 207 . Google Scholar CrossRef Search ADS PubMed Laanisto , L. and Niinemets , Ü. 2015 Polytolerance to abiotic stresses: how universal is the shade-drought tolerance trade-off in woody species . Glob. Ecol. Biogeogr. 24 , 571 – 580 . Google Scholar CrossRef Search ADS PubMed Langsaeter , A. 1941 Om tynning i enaldret gran-og furuskog . Meddel. f. d. Nor. Skogforsoksves 8 , 131 – 216 . Lawson , J.R. , Fryirs , K.A. and Leishman , M.R. 2015 Hydrological conditions explain variation in wood density in riparian plants of south-eastern Australia . J. Ecol. 103 , 945 – 956 . doi:10.1111/1365-2745.12408 . Google Scholar CrossRef Search ADS Long , J.N. 1985 A practical approach to density management . For. Chr. 61 , 23 – 27 . Google Scholar CrossRef Search ADS Long , J.N. and Daniel , T.W. 1990 Assessment of growing stock in uneven-aged stands . West. J. Appl. For. 5 , 93 – 96 . Long , J.N. , Dean , T.J. and Roberts , S.D. 2004 Linkages between silviculture and ecology: examination of several important conceptual models . For. Eco. Manag. 200 , 249 – 261 . Google Scholar CrossRef Search ADS Long , J.N. and Smith , F.W. 1984 Relation between size and density in developing stands: a description and possible mechanisms . For. Ecol. Manag. 7 , 191 – 206 . doi:10.1016/0378-1127(84)90067-7 . Google Scholar CrossRef Search ADS Long , J.N. and Vacchiano , G. 2014 A comprehensive framework of forest stand property-density relationships: perspectives for plant population ecology and forest management . Ann. For. Sci. 71 , 325 – 335 . doi:10.1007/s13595-013-0351-3 . Google Scholar CrossRef Search ADS Lonsdale , W.M. 1990 The self-thinning rule: dead of alive? Ecology 71 , 1373 – 1388 . Google Scholar CrossRef Search ADS Niinemets , Ü. and Valladares , F. 2006 Tolerance to shade, drought, and waterlogging of temperate northern hemisphere trees and shrubs . Ecol. Monogr. 76 , 521 – 547 . doi:10.1890/0012-9615(2006)076[0521:TTSDAW]2.0.CO;2 . Google Scholar CrossRef Search ADS Oliver , C.D. and Larson , B.C. 1996 Forest Stand Dynamics . John Wiley and Sons . Pompa-García , M. and Venegas-González , A. 2016 Temporal variation of wood density and carbon in two elevational sites of Pinus cooperi in relation to climate response in northern Mexico . PLoS One 11 , e0156782 . doi:10.1371/journal.pone.0156782 . Google Scholar CrossRef Search ADS PubMed Pretzsch , H. 2014 Canopy space filling and tree crown morphology in mixed-species stands compared with monocultures . For. Ecol. Manag. 327 , 251 – 264 . doi:10.1016/j.foreco.2014.04.027 . Google Scholar CrossRef Search ADS Pretzsch , H. and Biber , P. 2005 A re-evaluation of Reineke’ s rule and stand density index . For. Sci. 51 , 304 – 320 . Pretzsch , H. and Biber , P. 2016 Tree species mixing can increase maximum stand density . Can. J. For. Res. 46 , 1179 – 1193 . doi:10.1139/cjfr-2015-0413 . Google Scholar CrossRef Search ADS Pretzsch , H. , del Río , M. , Ammer , C. , Avdagic , A. , Barbeito , I. , Bielak , K. , et al. 2015 Growth and yield of mixed versus pure stands of Scots pine (Pinus sylvestris L.) and European beech (Fagus sylvatica L.) analysed along a productivity gradient through Europe . Eur. J. For. Res. 134 , 927 – 947 . doi:10.1007/s10342-015-0900-4 . Google Scholar CrossRef Search ADS Pretzsch , H. and Mette , T. 2008 Linking stand-level self-thinning allometry to the tree-level leaf biomass allometry . Trees 22 , 611 – 622 . doi:10.1007/s00468-008-0231-x . Google Scholar CrossRef Search ADS Pretzsch , H. and Schütze , G. 2015 Effect of tree species mixing on the size structure, density, and yield of forest stands . Eur. J. For. Res. 135 , 1 – 22 . doi:10.1007/s10342-015-0913-z . Google Scholar CrossRef Search ADS R-Core Team , 2016 . R: a language and environment for statistical computing. Reineke , L.H. 1933 Perfecting a stand-density-index for even-age forests . J. Agric. Res. 46 , 627 – 638 . Reyes-Hernandez , V. , Comeau , P.G. and Bokalo , M. 2013 Static and dynamic maximum size–density relationships for mixed trembling aspen and white spruce stands in western Canada . For. Ecol. Manag. 289 , 300 – 311 . doi:10.1016/j.foreco.2012.09.042 . Google Scholar CrossRef Search ADS Valbuena , P. , Peso , C. and Bravo , F. 2008 Stand density management diagrams for two mediterranean pine species in Eastern Spain . Invest. Agrar. Sist. Rec. 17 , 97 – 104 . Google Scholar CrossRef Search ADS Valladares , F. and Niinemets , U. 2008 Shade tolerance, a key plant feature of complex nature and consequences . Ann. Rev. Ecol. Evol. Syst. 39 , 237 – 257 . Google Scholar CrossRef Search ADS Vincent , G. and Harja , D. 2008 Exploring ecological significance of tree crown plasticity through three-dimensional modelling . Ann. Bot. 101 , 1221 – 1231 . doi:10.1093/aob/mcm189 . Google Scholar CrossRef Search ADS PubMed Vospernik , S. and Sterba , H. 2014 Do competition-density rule and self-thinning rule agree? Ann. For. Sci. 72 , 379 – 390 . doi:10.1007/s13595-014-0433-x . Google Scholar CrossRef Search ADS Weiskittel , A. , Gould , P. and Temesgen , H. 2009 Sources of variation in the self-thinning boundary line for three species with varying levels of shade tolerance . For. Sci. 55 , 84 – 93 . Woodall , C.W. , Miles , P.D. and Vissage , J.S. 2005 Determining maximum stand density index in mixed species stands for strategic-scale stocking assessments . For. Ecol. Manag. 216 , 367 – 377 . doi:10.1016/j.foreco.2005.05.050 . Google Scholar CrossRef Search ADS Wright , I.J. , Ackerly , D.D. , Bongers , F. , Harms , K.E. , Ibarra-Manriquez , G. , Martinez-Ramos , M. , et al. 2007 Relationships among ecologically important dimensions of plant trait variation in seven neotropical forests . Ann. Bot. 99 , 1003 – 1015 . doi:10.1093/aob/mcl066 . Google Scholar CrossRef Search ADS PubMed Zeide , B. 1987 Analysis of the 3/2 power law of self-thinning . For. Sci. 33 , 517 – 537 . Zeide , B. 2004 Optimal stand density: a solution . Can. J. For. Res. 34 , 846 – 854 . doi:10.1139/X03-258 . Google Scholar CrossRef Search ADS Zeller , L. , Ammer , C. , Annighöfer , P. , Biber , P. , Marshall , J. , Schütze , G. , et al. 2017 Tree ring wood density of Scots pine and European beech lower in mixed-species stands compared with monocultures . For. Ecol. Manag. 400 , 363 – 374 . doi:10.1016/j.foreco.2017.06.018 . Google Scholar CrossRef Search ADS Zhang , L. , Bi , H. , Gove , J.H. and Heath , L.S. 2005 A comparison of alternative methods for estimating the self-thinning boundary line . Can. J. For. Res. 35 , 1507 – 1514 . doi:10.1139/x05-070 . Google Scholar CrossRef Search ADS Zhang , J. , Oliver , W.W. and Powers , R.F. 2013 Reevaluating the self-thinning boundary line for ponderosa pine (Pinus ponderosa) forests . Can. J. For. Res. 43 , 971 . doi:10.1139/cjfr-2013-0133 . © Institute of Chartered Foresters, 2018. All rights reserved. For Permissions, please e-mail: journals.permissions@oup.com. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Forestry: An International Journal Of Forest Research Oxford University Press

Maximum stand density strongly depends on species-specific wood stability, shade and drought tolerance

Loading next page...
 
/lp/ou_press/maximum-stand-density-strongly-depends-on-species-specific-wood-BBEZ0i8kxl
Publisher
Institute of Chartered Foresters
Copyright
© Institute of Chartered Foresters, 2018. All rights reserved. For Permissions, please e-mail: journals.permissions@oup.com.
ISSN
0015-752X
eISSN
1464-3626
D.O.I.
10.1093/forestry/cpy006
Publisher site
See Article on Publisher Site

Abstract

Abstract The abundance of stems in crowded populations and the subsequent self-thinning is a key issue in forest stand dynamics. However, the mechanisms that control self-thinning are challenging to model. Although some attempts to include climate and structural traits like specific gravity (SG) are promising, they remain confined to North American species. In this study we aimed to disentangle how SG along with two major abiotic stress tolerances, i.e. shade and drought tolerance, contribute to the maximum density of a forest stand across a climatic gradient in Europe, and thus test the validity of the species-specific trait control over stand density. We propose a modelling approach that incorporates the tolerance to drought and shade in the determination of maximum relative stocking. Here, relative stocking refers to the degree of tree crowding in forest ecosystems. A relative stocking base model where specific wood density is inversely related to stand density is modified, adding normalized indices of drought and shade stress tolerance. We used available species tolerance rankings modulated by stress intensity to analyse the effects of abiotic stress polytolerance or trade-offs in the study area which represent an environmental gradient from Alpine to Mediterranean climate in northern Spain. Results indicated that the role of drought tolerance in controlling maximum stand density is stronger in warmer sites. The simultaneous tolerance to shade and drought results in less carrying capacity of sites. In those sites where there is no water limitation but minimum temperature is very low the tolerance to bending stress (i.e. specific gravity) explains better the maximum tree occupancy. Introduction Forest stand density is a key variable in forestry and ecology. Understanding the reasons why tree species occupy an area and the maximum capacity of each site to harbour a certain number of individuals is at the core of management decisions. Early silvicultural studies highlighted the concept that stand density and its control allows managers to keep production within optimum levels to compensate the growth trade-off between stand- and tree-level (Langsaeter, 1941; Assmann, 1971). Long (1985) summarized this hypothesis establishing the density limits between the thresholds of lower full site occupancy and the onset of self-thinning. These limits are based on the maximum stand density of a site that represents the maximum degree of site occupancy (Dean and Baldwin, 1996). The determination of maximum density has largely been explored in terms of maximum size-density relationships (MSDR). Most forestry studies explore MSDR using Reineke’s equation (Reineke, 1933) and the associated Stand Density Index (SDI) that relates the number of stems as a function of the size of the trees in relation to a fixed diameter, usually 25 cm (Pretzsch and Biber, 2005). Similarly, ecological studies have dealt with the −3/2 self-thinning law postulated by Yoda et al. (1963) that relates the individual mean volume with the number of individuals. Both relationships in their original formulation are based on the assumption that there is a limit in the relationship between size and density and they are equivalent if tree volume is proportional to quadratic mean diameter (Burkhart and Tomé, 2012). Maximum forest stand density represents an upper limit to the occupancy of a site, and growth is only possible after the death of some individuals. This upper limit on potential site occupancy has been considered to be species- and site-specific by several authors using different variables to characterize the stand. Condés et al. (2017) used de Martonne’s aridity index (de Martonne, 1926) for Scots pine (Pinus sylvestris L.) and European beech (Fagus sylvatica L.) whereas Zhang et al. (2013) included site index to model MSDR for ponderosa pine (Pinus ponderosa). The models are fitted for specific species compositions and the results cannot be applied to different compositions. Thus, while potentially useful for well-defined management situations, the ecological inferences from such models remain rather limited. The simultaneous representation of density boundary lines across many different species, growing in both pure and mixed stands and in different environments is still seldom explored. Site index was also used in conjunction with other stand-specific attributes like site index, climate, stand purity and origin to fit a MSDR for different forest populations and communities of Douglas-fir (Pseudotsuga menziesii var. menziesii [Mirb.] Franco), western hemlock (Tsuga heterophylla (Raf.) Sarg.) and red alder (Alnus rubra Bong.) (Weiskittel et al. 2009). Similarly, site quality and species composition were included in modelling maximum density of trembling aspen (Populus tremuloides Michx.) and white spruce (Picea glauca (Moench) Voss.) stands (Reyes-Hernandez et al., 2013). Another approach is that based on species-specific traits like specific wood gravity, SG (Woodall et al., 2005; Ducey and Knapp, 2010) or a combination of two traits, like SG and shade tolerance (ST) modulated by climate (Ducey et al., 2017). The hypothesis behind the use of species-specific traits to predict maximum density is based on stress tolerance. Specific wood density is associated with tolerance to bending stress, as developed by Dean and Baldwin (1996) who first postulated that species with light wood would need more stems to support the same foliage mass per unit area than species with harder wood. The inverse relationship between maximum stand density, defined as Reineke’s stand density and SG was presented by Woodall et al. (2005) whose relationship can be applied to mixtures as long as specific gravity is known for all component species. Ducey and Knapp (2010) combined Woodall et al.’s (2005) approach with the summation method to apply the additive stand density index (ASDI) developed by Long and Daniel (1990) in heterogeneous forests along with the Curtis’ (1971) modified version of the tree area ratio (Chisman and Schumacher, 1940). The resulting model predicts the relative stocking of a stand based on the specific gravity of tree species. Recently, Ducey et al. (2017) added tolerance to low light environments in addition to bending stress. Their results showed how climate variables modulated stress tolerance and provided the improvement to a more mechanistic explanation of site occupancy. Shade tolerance is a key mechanism shaping forest communities (Valladares and Niinemets, 2008) and it has been used to estimate maximum stand density (Ducey et al., 2017). Other abiotic stresses might be more important, however, in other climatic regions like the Mediterranean where drought events are more harmful than low light levels. Thus, we must test the idea that a single or two-trait approach may be insufficient to characterize maximum density and raise the issue of dealing with polytolerance, i.e. the ability of a species to simultaneously tolerate different abiotic stressors. In this paper we aimed to identify which trait contributes more to the maximum density of a forest stand, and if the relative contribution of two major types of abiotic stress tolerance, i.e. shade and drought tolerance, changes across a climatic gradient. Our working hypotheses are: Mechanical properties, like the wood specific gravity, have an inverse relationship with stocking of forest stands across a climatic gradient. Shade tolerance value, as modulated by cold temperature, affects maximum density inversely: the more shade tolerant the species and colder the site the less number of individuals needed to occupy the stand. Drought tolerance is more important in drought-prone areas, like in Mediterranean areas, where fewer individuals are needed to fully occupy the stand due to site conditions. Methods Study area, stand and climate data The study is carried out in a transition zone covering Mediterranean, Atlantic, and Alpine climate in Navarra province, northern Spain (Figure 1). The mean annual temperature is 10.5°C, ranging from 4.3 to 14.2°C, and annual rainfall ranges from 373 to 2629 mm, averaging 1378 mm (Gonzalo-Jiménez, 2010). Stand information comprising species composition, number of individuals by unit area (N, stems ha−1), and quadratic mean diameter (dg, cm) were extracted from the Third Spanish National Forest Inventory (NFI) available at www.mapama.org. Each NFI plot consists of four concentric subplots with radius 5, 10, 15 and 25 m where trees with DBH equal to or over 7.5, 12.5, 22.5 and 42.5 cm were measured, respectively. To get the number of stems by unit area an expansion factor is applied for each subplot, i.e. for the 10 m subplot each tree represents 10 000/(π × 102) individuals. We removed from the database plots in fertilized or irrigated Populus x canadensis plantations to make sure that intensive management did not affect the intensity of abiotic stress. We retained 2803 plots in four climatic regions built upon Allué-Andrade’s (1990) phytoclimatic classification: Mediterranean, Mountainous, Maritime and Transition (see Supplementary Table 1 for description of climatic regions). A species-specific wood gravity for each species were retrieved from Gutiérrez and Plaza (1967) whereas shade and drought tolerance (ST and DT, respectively) was obtained from Niinemets and Valladares (2006). For those species without a tolerance value we applied the average value for the genus. Table 1 shows main stand level properties and stress tolerance range by climatic gradient. The dataset comprises 66 species ranging from temperate to Mediterranean-type forests. The most abundant species are European beech (Fagus sylvatica L), Scots pine (Pinus sylvestris L.), Holm oak (Quercus ilex L.), oak (Q. robur and Q. petraea) and quejigo oak (Quercus faginea Lam.) (Supplementary Table 2 describes the main dendrometric characteristic of each species growing in the study region). Figure 1 View largeDownload slide Study area (Navarra, Spain), biogeoclimatic zones and National Forest Inventory plots used in this study. Figure 1 View largeDownload slide Study area (Navarra, Spain), biogeoclimatic zones and National Forest Inventory plots used in this study. Table 1 Descriptive stand variables, specific gravity and stress tolerances by bioregions in the study area. Bioregion # plots Hdom (m) BA (m2 ha−1) dg (cm) N (stems ha−1) SG ST DT Maritime 1044 17 (5.3) 3.5–33 28.1 (10.6) 0.7–78.2 26.1 (11.5) 8–111.5 755.6 (577) 5.1–3738.4 0.7 (0.1) 0.4–1 3.4 (1.1) 1.2–4.6 3 (0.8) 1.2–5 Mediterranean 374 12.4 (6) 3–27.6 22.8 (11.5) 0.4–52.6 19 (7) 7.6–81.3 998.8 (677.9) 10.2–2910.8 0.7 (0.2) 0.4–1 2.4 (0.8) 1.4–4.6 3.6 (1.1) 1.2–5 Mountainous 1121 17.8 (6) 2.5–36.5 31.5 (13.8) 0.4–82.7 25.6 (10) 7.6–95.5 808.2 (561.8) 5.1–3243.2 0.7 (0.1) 0.4–1 3.2 (1.4) 1.2–4.6 3.3 (0.9) 1.7–5 Transition 264 10 (4.1) 3–28.3 22.6 (10.4) 0.6–56.1 17.7 (8.7) 7.7–113.9 1214.2 (775.2) 10.2–3448.4 0.9 (0.2) 0.4–1 2.9 (0.8) 1.4–4.6 3.9 (0.9) 1.7–5 Bioregion # plots Hdom (m) BA (m2 ha−1) dg (cm) N (stems ha−1) SG ST DT Maritime 1044 17 (5.3) 3.5–33 28.1 (10.6) 0.7–78.2 26.1 (11.5) 8–111.5 755.6 (577) 5.1–3738.4 0.7 (0.1) 0.4–1 3.4 (1.1) 1.2–4.6 3 (0.8) 1.2–5 Mediterranean 374 12.4 (6) 3–27.6 22.8 (11.5) 0.4–52.6 19 (7) 7.6–81.3 998.8 (677.9) 10.2–2910.8 0.7 (0.2) 0.4–1 2.4 (0.8) 1.4–4.6 3.6 (1.1) 1.2–5 Mountainous 1121 17.8 (6) 2.5–36.5 31.5 (13.8) 0.4–82.7 25.6 (10) 7.6–95.5 808.2 (561.8) 5.1–3243.2 0.7 (0.1) 0.4–1 3.2 (1.4) 1.2–4.6 3.3 (0.9) 1.7–5 Transition 264 10 (4.1) 3–28.3 22.6 (10.4) 0.6–56.1 17.7 (8.7) 7.7–113.9 1214.2 (775.2) 10.2–3448.4 0.9 (0.2) 0.4–1 2.9 (0.8) 1.4–4.6 3.9 (0.9) 1.7–5 Average, standard deviation in parenthesis and range, minimum-maximum for Hdom: Dominant height, BA: basal area, dg: quadratic mean diameter, N: number of stems per hectare, SG: specific gravity according to Gutierrez Oliva and Plaza Pulgar (1967), ST: averaged shade tolerance, and DT: drought tolerance according to Niinimets and Valladares (2008). The dominant species and its respective SG, ST and DT values by bioregions are European beech (0.75, 4.56, 2.4) and Scots pine (0.50, 1.67, 4.34) in Maritime and Mountainous regions, Black pine (0.58, 2.1, 4.38) and Aleppo pine (0.55, 1.35, 4.97) in the Mediterranean region and Holm oak (1.00, 3.02, 4.72) and Faginea oak (0.95, 2.73, 3.02) in the Transition bioregion. Table 1 Descriptive stand variables, specific gravity and stress tolerances by bioregions in the study area. Bioregion # plots Hdom (m) BA (m2 ha−1) dg (cm) N (stems ha−1) SG ST DT Maritime 1044 17 (5.3) 3.5–33 28.1 (10.6) 0.7–78.2 26.1 (11.5) 8–111.5 755.6 (577) 5.1–3738.4 0.7 (0.1) 0.4–1 3.4 (1.1) 1.2–4.6 3 (0.8) 1.2–5 Mediterranean 374 12.4 (6) 3–27.6 22.8 (11.5) 0.4–52.6 19 (7) 7.6–81.3 998.8 (677.9) 10.2–2910.8 0.7 (0.2) 0.4–1 2.4 (0.8) 1.4–4.6 3.6 (1.1) 1.2–5 Mountainous 1121 17.8 (6) 2.5–36.5 31.5 (13.8) 0.4–82.7 25.6 (10) 7.6–95.5 808.2 (561.8) 5.1–3243.2 0.7 (0.1) 0.4–1 3.2 (1.4) 1.2–4.6 3.3 (0.9) 1.7–5 Transition 264 10 (4.1) 3–28.3 22.6 (10.4) 0.6–56.1 17.7 (8.7) 7.7–113.9 1214.2 (775.2) 10.2–3448.4 0.9 (0.2) 0.4–1 2.9 (0.8) 1.4–4.6 3.9 (0.9) 1.7–5 Bioregion # plots Hdom (m) BA (m2 ha−1) dg (cm) N (stems ha−1) SG ST DT Maritime 1044 17 (5.3) 3.5–33 28.1 (10.6) 0.7–78.2 26.1 (11.5) 8–111.5 755.6 (577) 5.1–3738.4 0.7 (0.1) 0.4–1 3.4 (1.1) 1.2–4.6 3 (0.8) 1.2–5 Mediterranean 374 12.4 (6) 3–27.6 22.8 (11.5) 0.4–52.6 19 (7) 7.6–81.3 998.8 (677.9) 10.2–2910.8 0.7 (0.2) 0.4–1 2.4 (0.8) 1.4–4.6 3.6 (1.1) 1.2–5 Mountainous 1121 17.8 (6) 2.5–36.5 31.5 (13.8) 0.4–82.7 25.6 (10) 7.6–95.5 808.2 (561.8) 5.1–3243.2 0.7 (0.1) 0.4–1 3.2 (1.4) 1.2–4.6 3.3 (0.9) 1.7–5 Transition 264 10 (4.1) 3–28.3 22.6 (10.4) 0.6–56.1 17.7 (8.7) 7.7–113.9 1214.2 (775.2) 10.2–3448.4 0.9 (0.2) 0.4–1 2.9 (0.8) 1.4–4.6 3.9 (0.9) 1.7–5 Average, standard deviation in parenthesis and range, minimum-maximum for Hdom: Dominant height, BA: basal area, dg: quadratic mean diameter, N: number of stems per hectare, SG: specific gravity according to Gutierrez Oliva and Plaza Pulgar (1967), ST: averaged shade tolerance, and DT: drought tolerance according to Niinimets and Valladares (2008). The dominant species and its respective SG, ST and DT values by bioregions are European beech (0.75, 4.56, 2.4) and Scots pine (0.50, 1.67, 4.34) in Maritime and Mountainous regions, Black pine (0.58, 2.1, 4.38) and Aleppo pine (0.55, 1.35, 4.97) in the Mediterranean region and Holm oak (1.00, 3.02, 4.72) and Faginea oak (0.95, 2.73, 3.02) in the Transition bioregion. Model structure The hypothesis was tested with a modified version of the model structure developed by Ducey et al. (2017). This model is based on the tree area ratio proposed by Chisman and Schumacher (1940) who defined the area occupied by the crown projection, i.e. the growing space requirement of a tree, as a quadratic function of stem diameter (equation (1)): Ai=c0+c1DBHi+c2DBHi2, (1) where Ai is the available growing space and DBHi is diameter at breast height of tree i, c0-c2 are parameters to be estimated. Summing all individual growing spaces results in the Tree Area Ratio (TAR, equation (2)) which represents the growing space occupied by trees in a stand. TAR=c0N+c1∑iDBHi+c2∑iDBHi2 (2) Although Chisman and Schumacher (1940) identified TAR with the potential crown space occupancy of a stand, other researchers have extended the interpretation to reflect a more abstract resource-based growing space requirement (sensuOliver and Larson, 1996). Setting TAR to 1 and estimating parameters c0,c1 and c2 with data from fully stocked stands, equation (2) can be used to predict the relative density of new stands. Ducey and Knapp (2010) took this approach but changed the growing space requirement to a function of the additive stand density index (Long and Daniel, 1990) combined with wood specific gravity following Woodall et al. (2005). More recently Ducey et al. (2017) modified the approach to accommodate covariates such as climate and functional traits like shade tolerance, e.g., Ai=(b0+∑j=1TbjTij)(DBHi25)1.6, (3) where Ai is the growing space of tree i, Tij is any plant trait of tree i (such as specific gravity or shade tolerance), DBHi is diameter at breast height of tree i and b0 through bj are parameters to be estimated. Summing over individual trees provides a relative density similar to TAR (equation (4)): RD=b0∑i=1n(DBHi25)1.6+∑j=1TbjTij∑i=1n(DBHi25)1.6. (4) Ducey et al. (2017) further extended the model to include climate variables as well as species functional traits. An appropriate fitting procedure is needed to estimate the parameters in equation (4). If, as in Chisman and Schumacher (1940), RD were set to equal 1 for all stands, and then least-squares regression were used to estimate the parameters, there would be an implicit assumption that the stands included in the data had an average stocking level corresponding to full stocking. Thus, it would be necessary to pre-screen the data to select only fully stocked stands. However, other fitting approaches do not need such pre-screening, as illustrated by Ducey and Knapp (2010) and Ducey et al. (2017) (see below in the Fitting procedure section). Stress tolerance and stress intensity Testing our hypothesis needs comparable coefficients and a reasonable way to modulate the stress tolerance as a function of the stress intensity. Our modification of Ducey et al.’s (2017) model includes a normalization of Niinemets and Valladares (2006) shade and drought tolerance rankings to make them comparable to the specific gravity values found in the dataset. We scaled DT and ST ranked values to a maximum of 1 by dividing the score value by the maximum value found in the dataset, renaming them as nDT and nST, respectively. In so doing, we can assess the relative contribution of each tolerance to the maximum carrying capacity of the stand as all three traits have a maximum value of 1 (e.g. maximum SG is 1 g cm−3 for holm oak, Quercus ilex L.). Instead of using continuous climatic variables to assess the climate effect on stress tolerance we built multipliers which reduced the normalized stress tolerance in harsh conditions, hereafter referred as stress tolerance modulators, allowing us to order species tolerance from 0 (no stress and the tolerance value does not influence the maximum density) to the maximum tolerance, i.e. the maximum normalized tolerance value nDT or nST. The contribution of DT to predict the maximum density is regulated by a multiplier that maximizes the tolerance from lower to higher drought stress levels: dI=kmax(k), (5) where dI is the drought intensity stress modulator of nDT, k is the drought intensity, measured for each plot as the area confined between the temperature and precipitation curves in a Walter-Lieth’s climodiagram, max(k) is the maximum drought intensity for the study area, this modulator ranges from 0, i.e. sites without drought period where the contribution of drought tolerance is negligible (k = 0), to 1, i.e. sites with maximum drought intensity (k = max(k)). Analogously, shade tolerance is modulated by minimum temperature, which is correlated with the length of the growing season (Laanisto and Niinemets, 2015; Ducey et al., 2017): cI=1−Tmin−minTminmaxTmin−minTmin, (6) where cI is the cold-stress modulator of nST, Tmin is the minimum temperature of the coldest month, minTmin and maxTmin are the minimum and the maximum of the minimum temperature of the coldest month, respectively, this modulator ranges from 1 (minimum minT, associated to cold sites where shade tolerance is maximum for the species) to 0 (warmer sites). Although specific gravity might vary, among others, with altitude and dry conditions (Chave et al., 2006; Pompa-García and Venegas-González, 2016) or exposure to stressors, like flooding (Lawson et al., 2015) we assume, according to our working hypothesis 1, that bending stress is always occurring across the species distribution as it determines the maximum amount of foliage that a stand can support (Dean and Baldwin, 1996). In addition, there is no significant relationship between specific gravity and site aridity measured by de Martonne Index in the study area (see Supplementary Figure 1). The trade-offs and synergies of stress tolerance modulators were analysed with dispersion plots of two contrasting broadleaves (Fagus sylvatica L. and Quercus ilex L.) and conifers (Pinus sylvestris L. and P. halepensis Mill.). A species with great variability in both multipliers, i.e. large dispersion in both axes, would indicate a large suitable habitat where there could be a shift in the relative importance of shade or drought tolerance and even polytolerance, whereas low variability in either of both multipliers would indicate absence of polytolerance and that stress tolerance to a single stressor is more important for the species to occupy the site. Hypothesis testing We started the fitting procedure from the base model with SG as a single variable as in Ducey and Knapp (2010), equation (4) (hypothesis 1), and adding one more normalized stress tolerance alternatively (nST and nDT) multiplied by its corresponding stress modulator (cI and dI, respectively). The subsequent models used for the rest of hypotheses are as follows: Testing hypothesis 2 expand equation (4) with a term including the tolerance to shade RD=a0∑iEFi(DBHi25)1.6+a1∑iEFiSGi(DBHi25)1.6+a2∑iEFi⋅(cI⋅nSTi)⋅(DBHi25)1.6, (7) where RD is relative density, EFi is the expansion factor applicable to each individual according to its position in the concentric plot design of the NFI, DBHi is the diameter of individual, SGi is the species-specific wood gravity for each individual, nSTi is the species-specific normalized shade tolerance for each individual, cI is the cold temperature stress modulator, and a0, a1 and a2 are parameters to be estimated. Hypothesis 3 was tested by adding a term for drought tolerance to equation (4): RD=a0∑iEFi(DBHi25)1.6+a1∑iEFiSGi(DBHi25)1.6+a3∑iEFi⋅(dI⋅nDTi)(DBHi25)1.6, (8) where Di is the drought stress modulator, nDTi is the species-specific normalized drought tolerance and the rest of symbols are as in equation (7). Alternatively, we tested the simultaneous relative importance of tolerance to all stressors: RD=a0∑iEFi(DBHi25)1.6+a1∑iEFiSGi(DBHi25)1.6+a2∑iEFi⋅(cI⋅nSTi)⋅(DBHi25)1.6+a3∑iEFi⋅(dI⋅nDTi)(DBHi25)1.6, (9) where symbols were previously defined. Model specification described in equations (5) and (6) does not allow for simultaneous drought and shade tolerance values in neither cold nor warm sites (i.e. polytolerance). We also fitted the complementary value in the cold index (1-cI) to allow for testing shade–drought tolerance in warmer sites (equation (10)): RD=a0∑iEFi(DBHi25)1.6+a1∑iEFiSGi(DBHi25)1.6+a2∑iEFi⋅((1−cI)⋅nSTi)⋅(DBHi25)1.6+a3∑iEFi⋅(dI⋅nDTi)(DBHi25)1.6, (10) where all symbols were previously defined. Fitting procedure Traditional fitting procedures of maximum stocking need selection of plots at full stocking levels. However, screening data from national forest inventory plots is not efficient and rather subjective as the stocking level for RD = 1 (full site occupancy) is not truly known and there is evidence that plots in the data set do not correspond to RD = 1 (Ducey and Knapp, 2010), especially in large areas where environmental conditions and species traits can vary affecting the degree of occupancy. Plantations can be an exception and consequently they were removed from the data set, representing less than 1 per cent of the original dataset. Fitting models derived from equation (4) by ordinary least squares would require an independent and normal error term, and the regression procedure would select coefficients such that the mean value of the plots included in the dataset corresponded to RD = 1. However, our purpose here is to find a combination of coefficients such that the calculated RD value would fall below RD = 1 for nearly all plots; in other words, we seek a maximum density relationship at RD = 1 that would fall above all observations if the observations were free of measurement errors and other anomalies. To avoid subjective pre-screening procedures when using National Forest Inventory datasets, a default value of RD = 1 to all plots is hypothesized, and the models (equations (7) through (10)) are fit by quantile regression instead of ordinary least squares. The hypothetical assignment of RD = 1 to all plots does not mean that the actual density of all plots is at the maximum. Instead, it is giving a reference value of RD that should be consistent across all plots (no matter what their tree list, species traits, and climate might be). In such case, quantile regression allows the actual density of individual plots to vary freely, but it does enforce a model fit so that the reference RD = 1 actually falls at an extreme quantile of the distribution of RD values for individual plots, rather than the mean. By incorporating both species traits and climate variables and predictors, quantile regression ensures that the reference RD = 1, taken as a conditional quantile, remains consistently near the upper end of the distribution. Quantile regression fits a linear model through any value in data as long as the modeller selects a quantile. A quantile is a proportion, q∈[0,1] splitting the data into two proportions where the q observations are below the regression line and the 1-q are above. Maximum density needs a high quantile in order to assure that plots are at maximum density when on the modelled line; occasional observations may fall above the line due to measurement errors or unusual spatial configurations of trees relative to the plot boundaries. We tested four q values leaving 0.5 per cent, 1 per cent, 2.5 per cent and 5 per cent of observations above the maximum density boundary, i.e. we fitted quantiles 0.995, 0.99, 0.975 and 0.95. By setting RD = 1 to all plots and using the quantile regression we are removing all subjectivity associated to plot selection. Quantile regression has been successfully applied to model maximum boundary lines (Zhang et al., 2005; Ducey and Knapp, 2010; Vospernik and Sterba, 2014; Ducey et al., 2017). Schwarz’s Bayesian Information Criteria (BIC) within each quantile was used to assess fitting performance: BIC=−2⋅ln(L)−k⋅ln(n), (11) where L is the log-likelihood of the data, k is the number of parameters and n the number of observations. The lower BIC the better the fit. The significance level for parameter estimation was set to 0.05. All models were fitted using quantreg package of R (Koenker, 2016; R Core Team, 2016). Evaluation of maximum SDI for pure and mixed stands Maximum density predictions were evaluated graphically across environmental stress gradients (cold temperature and drought) between and within functional groups and species composition (broadleaves vs conifers and pure vs mixed-species stands). The maximum stand density relationship for two important commercial species, Scots pine (Pinus sylvestris L.) and European beech (Fagus sylvatica L.), has been previously assessed in the study area by Condés et al. (2017, 2013). Their maximum SDI was used as benchmark against to discuss our model’s output. In order to contrast results against climatic regions we also included comparisons of two Mediterranean species Aleppo pine (Pinus halepensis Mill.) and Holm oak (Quercus ilex L.) Additionally, we compared the SDI predicted by our model with information of SDI retrieved in the literature for species also found in our database: Pedunculate and Sessile oak (Quercus robur and Quercus petraea), black pine (Pinus nigra Arn.) and Silver fir (Abies alba L.) using values published by Pretzsch and Biber (2005) and Vospernik and Sterba (2014). Results Relationship and variability between normalized stress tolerance values Simultaneous variation in tolerance to drought and cold would indicate that a species or group of species can be found in a wider environmental range. Differences in polytolerance to shade and drought were detected by functional groups (Figure 2a and b). Plots of broadleaves growing in pure and mixed stands showed the highest variability for simultaneous normalized drought and shade tolerance whereas conifers did not show polytolerance either in pure or in mixed stands. Figure 2 View largeDownload slide Relationship between drought and shade tolerance ranking modulated by drought and cold temperature stress intensity (dI and cI). Large variability in both axes would indicate polytolerance. Figure 2 View largeDownload slide Relationship between drought and shade tolerance ranking modulated by drought and cold temperature stress intensity (dI and cI). Large variability in both axes would indicate polytolerance. A graphical analysis of stress modulators for two contrasting species per functional group showed that Scots pine and European beech are located mostly on drought free sites but with a large variability of minimum temperatures indicating that both species are not drought tolerant, whereas Holm oak (Quercus ilex L.) and Aleppo pine (Pinus halepensis Mill.) can be found over both a large range of drought intensity values and minimum temperature range (Figure 3a and b). Figure 3 View largeDownload slide (a) Relationship between cold and drought multipliers for broadleaves in Navarra (Spain). Large point dispersion in both axes would indicate simultaneous tolerance for cold and drought stress. Filled circles for Fagus sylvatica, void circles for Quercus ilex. (b) Relationship between cold and drought multipliers for conifers in Navarra (Spain). Large point dispersion in both axes would indicate simultaneous tolerance for cold and drought stress. Filled circles for Pinus sylvestris, void circles for Pinus halepensis. Figure 3 View largeDownload slide (a) Relationship between cold and drought multipliers for broadleaves in Navarra (Spain). Large point dispersion in both axes would indicate simultaneous tolerance for cold and drought stress. Filled circles for Fagus sylvatica, void circles for Quercus ilex. (b) Relationship between cold and drought multipliers for conifers in Navarra (Spain). Large point dispersion in both axes would indicate simultaneous tolerance for cold and drought stress. Filled circles for Pinus sylvestris, void circles for Pinus halepensis. Best model selection SG as the only explanatory variable (Table 2) was significant at all quantiles (equation (4), hypothesis 1 accepted) leading to an implied maximum density for Scots pine ranging from 1427 stems ha−1 at quantile 0.995 to 1121 stems ha−1 at quantile 0.95. European beech showed a maximum density between 1233 stems ha−1 and 945 stems ha−1 for the same quantiles, respectively. The addition of shade tolerance alone led to non-significant parameters at 0.05 significant level (equation (7), hypothesis 2 rejected), whereas drought tolerance alone led to significant models at all quantiles (equation (8), hypothesis 3 accepted). The full model with ST acting in cold sites and DT in warm sites (equation (9)) led to non-significant parameters for ST at all quantiles. The best fitted model in terms of BIC includes both normalized stress tolerance values modulated by temperature and drought intensity in warmer sites (equation (10)). This model at quantile 0.975 will be used for further SDI predictions as it represents the best fitted upper boundary line. Table 2 Estimates of relative density models including functional traits (SG, specific gravity; DT, drought tolerance; ST, shade tolerance) and abiotic stress modulators (iC cold temperature modulator and iD drought modulator). Model Variable q0.995 q0.99 q0.975 q0.95 Estimate (s.e.) BIC Estimate (s.e.) BIC Estimate (s.e.) BIC Estimate (s.e.) BIC (4) Intercept 0.00048 (0.00009) 3364.3 0.00043 (0.0001) 3305.9 0.00044 (0.00011) 3206.57 0.00051 (0.0001) 3113.18 SG 0.00044 (0.00013) 0.00059 (0.00016) 0.00068 (0.00015) 0.00073 (0.00013) (7) Intercept 0.00048 (0.00009) 0.0004 (0.0001) 0.00044 (0.00012) 0.00051 (0.0001) SG 0.00044 (0.00013) 0.00061 (0.00017) 0.00068 (0.00015) 0.00078 (0.00014) iC x ST n.s. n.s. n.s. n.s. (8) Intercept 0.00048 (0.00008) 3317.44 0.00043 (0.0001) 3273.17 0.00044 (0.00011) 3171.97 0.00051 (0.00009) 3077.7 SG 0.00044 (0.00013) 0.00059 (0.00016) 0.00068 (0.00016) 0.00071 (0.00013) iD x DT 0.00129 (0.00033) 0.00126 (0.00022) 0.0012 (0.00022) 0.00121 (0.00028) (9) Intercept 0.00048 (0.00008) 0.0004 (0.0001) 0.00014 (0.00012) 0.00051 (0.00009) SG 0.0044 (0.00011) 0.00061 (0.00018) 0.00068 (0.00015) 0.00076 (0.00013) iC x ST n.s. n.s. n.s. n.s. iD x DT 0.00129 (0.00039) 0.00127 (0.00023) 0.0012 (0.00025) 0.00115 (0.00026) (10) Intercept 0.00049 (0.00008) 0.00051 (0.00011) 0.00052 (0.0001) 3116.91 0.00054 (0.00009) 3035.1 SG n.s. n.s. 0.00033 (0.00016) 0.00045 (0.00014) (1−iC) x ST 0.00029 (0.00013) n.s. 0.00046 (0.00015) 0.00041 (0.00011) iD x DT 0.0013 (0.00021) 0.00128 (0.00023) 0.00122 (0.0002) 0.00115 (0.00029) Model Variable q0.995 q0.99 q0.975 q0.95 Estimate (s.e.) BIC Estimate (s.e.) BIC Estimate (s.e.) BIC Estimate (s.e.) BIC (4) Intercept 0.00048 (0.00009) 3364.3 0.00043 (0.0001) 3305.9 0.00044 (0.00011) 3206.57 0.00051 (0.0001) 3113.18 SG 0.00044 (0.00013) 0.00059 (0.00016) 0.00068 (0.00015) 0.00073 (0.00013) (7) Intercept 0.00048 (0.00009) 0.0004 (0.0001) 0.00044 (0.00012) 0.00051 (0.0001) SG 0.00044 (0.00013) 0.00061 (0.00017) 0.00068 (0.00015) 0.00078 (0.00014) iC x ST n.s. n.s. n.s. n.s. (8) Intercept 0.00048 (0.00008) 3317.44 0.00043 (0.0001) 3273.17 0.00044 (0.00011) 3171.97 0.00051 (0.00009) 3077.7 SG 0.00044 (0.00013) 0.00059 (0.00016) 0.00068 (0.00016) 0.00071 (0.00013) iD x DT 0.00129 (0.00033) 0.00126 (0.00022) 0.0012 (0.00022) 0.00121 (0.00028) (9) Intercept 0.00048 (0.00008) 0.0004 (0.0001) 0.00014 (0.00012) 0.00051 (0.00009) SG 0.0044 (0.00011) 0.00061 (0.00018) 0.00068 (0.00015) 0.00076 (0.00013) iC x ST n.s. n.s. n.s. n.s. iD x DT 0.00129 (0.00039) 0.00127 (0.00023) 0.0012 (0.00025) 0.00115 (0.00026) (10) Intercept 0.00049 (0.00008) 0.00051 (0.00011) 0.00052 (0.0001) 3116.91 0.00054 (0.00009) 3035.1 SG n.s. n.s. 0.00033 (0.00016) 0.00045 (0.00014) (1−iC) x ST 0.00029 (0.00013) n.s. 0.00046 (0.00015) 0.00041 (0.00011) iD x DT 0.0013 (0.00021) 0.00128 (0.00023) 0.00122 (0.0002) 0.00115 (0.00029) Table 2 Estimates of relative density models including functional traits (SG, specific gravity; DT, drought tolerance; ST, shade tolerance) and abiotic stress modulators (iC cold temperature modulator and iD drought modulator). Model Variable q0.995 q0.99 q0.975 q0.95 Estimate (s.e.) BIC Estimate (s.e.) BIC Estimate (s.e.) BIC Estimate (s.e.) BIC (4) Intercept 0.00048 (0.00009) 3364.3 0.00043 (0.0001) 3305.9 0.00044 (0.00011) 3206.57 0.00051 (0.0001) 3113.18 SG 0.00044 (0.00013) 0.00059 (0.00016) 0.00068 (0.00015) 0.00073 (0.00013) (7) Intercept 0.00048 (0.00009) 0.0004 (0.0001) 0.00044 (0.00012) 0.00051 (0.0001) SG 0.00044 (0.00013) 0.00061 (0.00017) 0.00068 (0.00015) 0.00078 (0.00014) iC x ST n.s. n.s. n.s. n.s. (8) Intercept 0.00048 (0.00008) 3317.44 0.00043 (0.0001) 3273.17 0.00044 (0.00011) 3171.97 0.00051 (0.00009) 3077.7 SG 0.00044 (0.00013) 0.00059 (0.00016) 0.00068 (0.00016) 0.00071 (0.00013) iD x DT 0.00129 (0.00033) 0.00126 (0.00022) 0.0012 (0.00022) 0.00121 (0.00028) (9) Intercept 0.00048 (0.00008) 0.0004 (0.0001) 0.00014 (0.00012) 0.00051 (0.00009) SG 0.0044 (0.00011) 0.00061 (0.00018) 0.00068 (0.00015) 0.00076 (0.00013) iC x ST n.s. n.s. n.s. n.s. iD x DT 0.00129 (0.00039) 0.00127 (0.00023) 0.0012 (0.00025) 0.00115 (0.00026) (10) Intercept 0.00049 (0.00008) 0.00051 (0.00011) 0.00052 (0.0001) 3116.91 0.00054 (0.00009) 3035.1 SG n.s. n.s. 0.00033 (0.00016) 0.00045 (0.00014) (1−iC) x ST 0.00029 (0.00013) n.s. 0.00046 (0.00015) 0.00041 (0.00011) iD x DT 0.0013 (0.00021) 0.00128 (0.00023) 0.00122 (0.0002) 0.00115 (0.00029) Model Variable q0.995 q0.99 q0.975 q0.95 Estimate (s.e.) BIC Estimate (s.e.) BIC Estimate (s.e.) BIC Estimate (s.e.) BIC (4) Intercept 0.00048 (0.00009) 3364.3 0.00043 (0.0001) 3305.9 0.00044 (0.00011) 3206.57 0.00051 (0.0001) 3113.18 SG 0.00044 (0.00013) 0.00059 (0.00016) 0.00068 (0.00015) 0.00073 (0.00013) (7) Intercept 0.00048 (0.00009) 0.0004 (0.0001) 0.00044 (0.00012) 0.00051 (0.0001) SG 0.00044 (0.00013) 0.00061 (0.00017) 0.00068 (0.00015) 0.00078 (0.00014) iC x ST n.s. n.s. n.s. n.s. (8) Intercept 0.00048 (0.00008) 3317.44 0.00043 (0.0001) 3273.17 0.00044 (0.00011) 3171.97 0.00051 (0.00009) 3077.7 SG 0.00044 (0.00013) 0.00059 (0.00016) 0.00068 (0.00016) 0.00071 (0.00013) iD x DT 0.00129 (0.00033) 0.00126 (0.00022) 0.0012 (0.00022) 0.00121 (0.00028) (9) Intercept 0.00048 (0.00008) 0.0004 (0.0001) 0.00014 (0.00012) 0.00051 (0.00009) SG 0.0044 (0.00011) 0.00061 (0.00018) 0.00068 (0.00015) 0.00076 (0.00013) iC x ST n.s. n.s. n.s. n.s. iD x DT 0.00129 (0.00039) 0.00127 (0.00023) 0.0012 (0.00025) 0.00115 (0.00026) (10) Intercept 0.00049 (0.00008) 0.00051 (0.00011) 0.00052 (0.0001) 3116.91 0.00054 (0.00009) 3035.1 SG n.s. n.s. 0.00033 (0.00016) 0.00045 (0.00014) (1−iC) x ST 0.00029 (0.00013) n.s. 0.00046 (0.00015) 0.00041 (0.00011) iD x DT 0.0013 (0.00021) 0.00128 (0.00023) 0.00122 (0.0002) 0.00115 (0.00029) According to the best model, the relative contribution of shade and drought tolerance to explain the maximum density for a given species increases in warm sites. Specific wood gravity and its relation to bending stress contributes more to explain maximum density in colder sites. Predicted maximum SDI values Differences were found in maximum density between phylogenetic groups. Gymnosperms showed higher maximum number of stems per unit area at a given diameter than angiosperms. Within each phylogenetic group drought tolerant species had lower density than drought intolerant species whereas shade intolerant species had higher density. The maximum density is reduced across the environmental stress gradient with the abiotic stress factor changing from the south where drought is more intense to north, where cold temperatures is the limiting factor. The maximum SDI for Scots pine ranges from 1244 to 1416 stems ha−1, whereas European beech maximum density ranged from 826 to 1191 stems ha−1 (Figure 4a). Drought tolerant species needed fewer stems to full occupy the stand. Aleppo pine ranged from 473 to 1225 stems ha−1 and Holm oak ranged from 725 to 946 stems ha−1 (Figure 4b). Figure 4 View largeDownload slide (a) Predicted maximum stand density index (stems ha−1) across temperature stress gradient in Navarra, Spain, for cold-tolerant/drought intolerant species. Grey circles for Pinus sylvestris and black circles for Fagus sylvatica. (b) Predicted maximum stand density index (stems ha−1) across drought stress gradient in Navarra, Spain, for cold-intolerant/drought tolerant species. Grey circles for Pinus halepensis and black circles for Quercus ilex. Figure 4 View largeDownload slide (a) Predicted maximum stand density index (stems ha−1) across temperature stress gradient in Navarra, Spain, for cold-tolerant/drought intolerant species. Grey circles for Pinus sylvestris and black circles for Fagus sylvatica. (b) Predicted maximum stand density index (stems ha−1) across drought stress gradient in Navarra, Spain, for cold-intolerant/drought tolerant species. Grey circles for Pinus halepensis and black circles for Quercus ilex. Discussion Maximum density is important to our understanding of species stand dynamics as it indicates the maximum occupancy or full-resource use by trees. This relationship has been described as the link between forestry and ecology (Jack and Long, 1996; Long et al., 2004) and it relies in the constancy of leaf area index (LAI) during self-thinning (Long and Smith, 1984). The assumption of constancy of LAI has helped to build the hypothesis that tolerance to bending stress is driven by structural traits like specific wood density. Woodall et al. (2005) demonstrated that under this assumption trees with lower wood density would need fewer individuals to support the same amount of foliage than species with higher wood density. This hypothesis could be behind the evidence that Scots pine and beech mixed stands in Europe have greater densities than pure stands (Pretzsch and Schütze, 2015; Pretzsch et al., 2015) as this type of mixture also has lower wood density (Zeller et al., 2017). So far, this hypothesis has been tested for North America species and we present here the first test of this hypothesis for European tree species expanding the validity of the stress-bending stress control over stand density. However, whether LAI remains constant during self-thinning has been put into question as foliar mass development is age-dependent indicating that the leaf mass that a stand can support is a function of stand development (Holdaway et al., 2008). These authors pointed that the change in leaf mass can be attributable to changes in crown vertical structure. On the other hand, Pretzsch and Mette (2008) found constant LAI in self-thinned stands of Norway spruce and beech but variable specific leaf area (SLA) with increasing diameter attributed to losses of shade leaves from suppressed dead trees. Self-thinning implies redistribution of the same amount of foliage in less and larger individuals (Long et al. 2004). However, the non-constancy of LAI is more evident in old stands where large aged trees cannot occupied the space available after mortality or thinning events leading to departures from a linear size-density trajectory (Long and Vacchiano, 2014). In general, the older the individual tree the less the foliage in high dense stands (Long et al. 2004). To reconcile the LAI constancy hypothesis with observations of age-dependent LAI, self-thinning boundary line modelling must combine the abiotic stress intensity and tolerance approach, as done in our model, and the existing relationship between crown plasticity and shade tolerance (Pretzsch, 2014; Jucker et al., 2015) or the competitive advantage of species with higher crown plasticity (Vincent and Harja, 2008). Equation (10) predicts a change in the relative contribution of stress tolerance from specific wood gravity in cold sites to increasing contribution of drought and shade tolerance in dry sites. This is in line with the finding of Condés et al. (2017) who suggested that plasticity in hydraulic structure outweighs mechanical control over stand density in arid environments. However, tree hydraulics might not be only relevant in dry areas as water transport has been also proposed to explain the inverse relationship between wood density and leaf mass in Neotropical forests (Wright et al., 2007). Abiotic stress tolerance is an important driver of woody plant distribution and productivity, as negative changes of environmental conditions may reduce the realized niche and potential biomass growth of a species (Laanisto and Niinemets, 2015). Polytolerance to stress is considered rare in nature leading to a trade-off between shade and drought tolerance (Niinemets and Valladares, 2006). Equation (9) in Table 2 would have yielded a lack of simultaneous effect of shade and drought tolerance in maximum density at both ends of a cold-warm spectrum (Supplementary Figure 2a). However, the lack of statistical significance for shade tolerance in equation (9) and the significance of its complementary value in equation (10) indicates that the relative importance of shade and drought tolerance increases in warmer sites (Supplementary Figure 2b). Tolerance to both stressors result in a reduction of the number of individuals needed to make full use of resources. However, species growing in warm sites are usually more drought tolerant than shade tolerant and the predicted reduction is mainly driven by drought stress, as it was reflected by the relative contribution of drought tolerance which is 2.7 higher than shade tolerance (Table 2). In contrast, species growing in cold sites showed an increasing effect of tolerance to bending stress on the number of individuals as compared to tolerance to shade and drought. Cold sites are expected to have a shorter growing season which has been pointed out to be a main reason for shade–drought trade-offs (Valladares and Niinemets, 2008; Laanisto and Niinemets, 2015). Equation (10) does not predict the presence or absence of polytolerance but it can be used to make inferences on maximum density about the effect of observed shade–drought trade-offs or polytolerance in a cold-warm spectrum. Warm sites would support less trees per unit area irrespective of tolerance ranking whereas a suite of polytolerant species would need fewer individuals than species with stronger trade-offs (Supplementary Figure 3). Linking this finding with the fact that polytolerance is associated with a higher probability of survival in harsh conditions, as suggested by Laanisto and Niinemets (2015), would have to be assessed in terms of constancy or stress tolerance dependence of the slope parameter in the density-size relationship. The basic model used in our approach keeps a constant slope in the density-size relationship (equation (4)). The constancy of slope in forestry and ecology has been a subject of intense debate (Lonsdale, 1990; Zeide, 2004, 1987) and findings suggest that the 1.6 constant slope is much more an exception than a rule (del Rı́o et al., 2001; Pretzsch and Biber, 2005). The mathematical computation of a variable slope according to environmental filters or stress tolerances in our model is challenging but it would improve our understanding of the mechanisms behind site occupancy. Nevertheless, our SDI predictions using stress intensity and tolerance are realistic and do not dramatically depart from predictions using a variable slope (del Rı́o et al., 2001; Pretzsch and Biber, 2005; Vospernik and Sterba, 2014; Condés et al., 2017). SDI for Scots pine in our study, 1307 stems ha−1, is consistent with 1297 stems ha−1 found in Condés et al. (2013) for the same region and slightly less than the 1442 stems ha−1 found by Río et al. (2008) for Spain. For European beech our model predicted less individuals, 966 stems ha−1, as compared with 1042 stems ha−1 found by Condés et al. (2013), nevertheless the range of our predictions embraced that of Condés’ as it ranged from 826 in colder sites to 1191 stems ha−1 in warmer sites. Our model predictions are closer to other estimations for beech in Europe (Figure 5), whereas SDI values for Silver fir (1134 stems ha−1) and Black pine (1246 stems ha−1) in Navarra are sightless less from those found in Austria, 1278 and 1424 stems ha−1, respectively (Vospernik and Sterba, 2014). Figure 5 View largeDownload slide Mean SDI (bars), standard error (whiskers), maximum SDI (triangles) found in the literature search for the species in the x-axis. Black points are predictions in Navarra using model 4 in Table 2. Figure 5 View largeDownload slide Mean SDI (bars), standard error (whiskers), maximum SDI (triangles) found in the literature search for the species in the x-axis. Black points are predictions in Navarra using model 4 in Table 2. Stress tolerance to shade is an important driver in shaping species distribution and community composition whereas structural traits, like wood density predict species’ competitive response (Kunstler et al., 2016). However, tolerance is not constant for a given species and it varies across genotypes and environmental conditions (Barnes et al., 1997). For example, shade tolerance is affected by nutrient and water availability (Valladares and Niinemets, 2008 and reference therein) and it decreases for Douglas-fir, western redcedar and hemlock with increasing soil moisture. Stress tolerance differs both with species ontogeny being more important in juvenile stages than in mature forests, and might differ during the same stage within a gradient of stress and species with a large distribution area. For example, Scots pine does not show the same shade tolerance across its distribution in Europe, being more shade intolerant in Central Europe than in the South. This variation in stress tolerance allows species to survive and grow in a wide range of light conditions. This gives a clue that scoring species according to stress tolerance is not an absolute value and the tolerance might be modulated by climatic factors as tested by Ducey et al. (2017) and normalized in this paper. Our model cannot be used to test the over- or under-density of mixed vs pure stands as the contribution of stress tolerance to maximum density is averaged across species, i.e. the resulting density cannot be higher that the individual contribution of each species. This is a caveat as long as several studies have shown higher density in mixed stands than in single-species ones (Pretzsch et al., 2015; Pretzsch and Biber, 2016) in association to changes in specific wood density in mixed stands (Zeller et al., 2017). A practical use of the model presented here is the substitution of stand density measures in empirical growth and yield model to include stress tolerance in the prediction of timber production. In so doing, the stand level management models will capture differences in yield production according to differences in stress tolerances. Silviculture aims to maintain density below the maximum density to keep stocking between the onset of competition and the onset of the density-dependent mortality where growth is constrained by resource competition and production loss due to mortality is not economically important. When management objectives promote individual tree growth the density is kept near the lower limit but if higher stand level yield is expected then the density is kept near the upper limit. These limits are set between 35 and 65 per cent of maximum density based on the seminal paper by Long (1985) and have been used in silvicultural prescriptions (Valbuena et al., 2008), however these values remain empirical and a mechanistic explanation is still lacking. Although our model does not provide mechanistic information about the exact limits it helps to accommodate the density taking into account the abiotic stress intensity and the species tolerance. For example, in the case of managed Scots pine stands in the study area the SDI limits should be set between 868 and 467 stems ha−1 in cold sites and between 544 and 293 stems ha−1 in the warmest sites. Conclusions We built a normalized approach to explain the maximum density needed to make full use of resources in terms of variability of the stress tolerance ranking according to Niinemets and Valladares (2006) modulated by the intensity of drought and low temperature stress. The model is based on the assumption that bending stress is concomitant to all sites in the study area assuming a constant LAI during self-thinning. A basic model of stocking degree inversely related to specific wood density is found significant in the study area showing an inverse relationship (hypothesis 1 accepted). The addition of normalized tolerance shade showed to be non-significant (hypothesis 2 rejected) whereas the combination of drought tolerance reduces the number of individuals that a site can support being this reduction stronger in warm sites (hypothesis 3 accepted). A species or group of species with high capacity to tolerate several abiotic stresses simultaneously, i.e. polytolerance, need far less individuals to fully occupy a site than species showing tolerance trade-offs. Linking polytolerance/trade-offs with a mechanistic explanation of survival at the stand level needs further investigation that should include variation due to stress tolerance in the slope parameter of the density-size relationship. Supplementary data Supplementary data are available at Forestry online. Acknowledgements The authors thank the anonymous reviewers and handling Editor for their constructive comments. Olivia L. Bartlett provided additional helpful suggestions on the manuscript. Funding The Spanish Minister of Economy and Competitiveness (grant number: AGL2014–51964-C2–2-R); the New Hampshire Agricultural Experiment Station (Scientific Contribution Number 2768); the USDA National Institute of Food and Agriculture McIntire-Stennis (Project Accession Number 1007007); and Short Term Scientific Mission provided by COST Action EuMIXFOR (COST Action number FP1206) to A.B.-O. Conflict of interest statement None declared. References Assmann , E . 1971 The Principles of Forest Yield Study . Pergamon press . Barnes , B.V. , Zak , D.R. , Denton , S.R. and Spurr , S.H. 1997 Forest Ecology . 4th ed . John Wiley and Sons . Burkhart , H.E. and Tomé , M. 2012 Modeling Forest Trees and Stands . Springer . Google Scholar CrossRef Search ADS Chave , J. , Muller-Landau , H.C. , Baker , T.R. , Easdlae , T.A. , Ter Steege , H. , et al. 2006 Regional and phylogenetic variation of wood denisty across 2456 neotropical tree species . Ecol. Appl. 16 , 2356 – 2367 . doi:10.3159/1095-5674(2007)134[301:SCROHD]2.0.CO;2 . Google Scholar CrossRef Search ADS PubMed Chisman , H.H. , and Schumacher , F.X. 1940 On the tree-area ratio and certain of its applications . J. For. 38 ( 4 ), 311 – 317 . Condés , S. , Del , M. and Sterba , H. 2013 Mixing effect on volume growth of Fagus sylvatica and Pinus sylvestris is modulated by stand density . For. Ecol. Manag. 292 , 86 – 95 . Google Scholar CrossRef Search ADS Condés , S. , Vallet , P. , Bielak , K. , Bravo-Oviedo , A. , Coll , L. , Ducey , M.J. , et al. 2017 Climate influences on the maximum size-density relationship in Scots pine (Pinus sylvestris L.) and European beech (Fagus sylvatica L.) stands . For. Ecol. Manag. 385 , 295 – 307 . doi:10.1016/j.foreco.2016.10.059 . Google Scholar CrossRef Search ADS Curtis , R.O. 1971 A tree area power function and related stand density measures for Douglas-fir . For. Sci. 17 , 146 – 159 . de Martonne , E. 1926 Une Nouvelle Fonction CIimatologique: L’Indice d’Aridite . La Météorol. 2 , 449 – 458 . Dean , T.J. and Baldwin , V.C. 1996 The relationship between Reineke’s stand-density index and physical stem mechanics . For. Ecol. Manag. 81 , 25 – 34 . doi:10.1016/0378-1127(95)03666-0 . Google Scholar CrossRef Search ADS del Rı́o , M. , Montero , G. and Bravo , F. 2001 Analysis of diameter–density relationships and self-thinning in non-thinned even-aged Scots pine stands . For. Ecol. Manag. 142 , 79 – 87 . doi:10.1016/S0378-1127(00)00341-8 . Google Scholar CrossRef Search ADS Ducey , M.J. and Knapp , R.A. 2010 A stand density index for complex mixed species forests in the northeastern United States . For. Ecol. Manag. 260 , 1613 – 1622 . doi:10.1016/j.foreco.2010.08.014 . Google Scholar CrossRef Search ADS Ducey , M.J. , Woodall , C.W. and Bravo-Oviedo , A. 2017 Climate and species functional traits influence maximum live tree stocking in the Lake States, USA . For. Ecol. Manag. 386 , 51 – 61 . doi:10.1016/j.foreco.2016.12.007 . Google Scholar CrossRef Search ADS Gonzalo-Jiménez , J. 2010 Diagnosis fitoclimática de la España Peninsular: hacia un modelo de clasificación funcional de la vegetación y de los ecosistemas peninsulares españoles. Edited by Organismo Autónomo de Parques Nacionales Gutiérrez , A. , and Plaza , F. 1967 Características fisico-mecánicas de las maderas Españolas. Madrid: Ministerio de Agricultura. Dirección General de Montes, Caza y Pesca Fluvial. IFIE . Holdaway , R.J. , Allen , R.B. , Clinton , P.W. , Davis , M.R. and Coomes , D.A. 2008 Intraspecific changes in forest canopy allometries during self-thinning . Funct. Ecol. 22 , 460 – 469 . doi:10.1111/j.1365-2435.2008.01388.x . Google Scholar CrossRef Search ADS Jack , S.B. and Long , J.N. 1996 Linkages between silviculture and ecology: an analysis of density management diagrams . For. Ecol. Manag. 86 , 205 – 220 . doi:10.1016/S0378-1127(96)03770-X . Google Scholar CrossRef Search ADS Jucker , T. , Bouriaud , O. and Coomes , D.A. 2015 Crown plasticity enables trees to optimize canopy packing in mixed-species forests . Funct. Ecol. 29 , 1078 – 1086 . doi:10.1111/1365-2435.12428 . Google Scholar CrossRef Search ADS Koenker , R. , 2016 . quantreg: Quantile Regression. R package version 5.29. Kunstler , G. , Falster , D. , Coomes , D.A. , Hui , F. , Kooyman , R.M. , Laughlin , D.C. , et al. 2016 Plant functional traits have globally consistent effects on competition . Nature 529 , 204 – 207 . Google Scholar CrossRef Search ADS PubMed Laanisto , L. and Niinemets , Ü. 2015 Polytolerance to abiotic stresses: how universal is the shade-drought tolerance trade-off in woody species . Glob. Ecol. Biogeogr. 24 , 571 – 580 . Google Scholar CrossRef Search ADS PubMed Langsaeter , A. 1941 Om tynning i enaldret gran-og furuskog . Meddel. f. d. Nor. Skogforsoksves 8 , 131 – 216 . Lawson , J.R. , Fryirs , K.A. and Leishman , M.R. 2015 Hydrological conditions explain variation in wood density in riparian plants of south-eastern Australia . J. Ecol. 103 , 945 – 956 . doi:10.1111/1365-2745.12408 . Google Scholar CrossRef Search ADS Long , J.N. 1985 A practical approach to density management . For. Chr. 61 , 23 – 27 . Google Scholar CrossRef Search ADS Long , J.N. and Daniel , T.W. 1990 Assessment of growing stock in uneven-aged stands . West. J. Appl. For. 5 , 93 – 96 . Long , J.N. , Dean , T.J. and Roberts , S.D. 2004 Linkages between silviculture and ecology: examination of several important conceptual models . For. Eco. Manag. 200 , 249 – 261 . Google Scholar CrossRef Search ADS Long , J.N. and Smith , F.W. 1984 Relation between size and density in developing stands: a description and possible mechanisms . For. Ecol. Manag. 7 , 191 – 206 . doi:10.1016/0378-1127(84)90067-7 . Google Scholar CrossRef Search ADS Long , J.N. and Vacchiano , G. 2014 A comprehensive framework of forest stand property-density relationships: perspectives for plant population ecology and forest management . Ann. For. Sci. 71 , 325 – 335 . doi:10.1007/s13595-013-0351-3 . Google Scholar CrossRef Search ADS Lonsdale , W.M. 1990 The self-thinning rule: dead of alive? Ecology 71 , 1373 – 1388 . Google Scholar CrossRef Search ADS Niinemets , Ü. and Valladares , F. 2006 Tolerance to shade, drought, and waterlogging of temperate northern hemisphere trees and shrubs . Ecol. Monogr. 76 , 521 – 547 . doi:10.1890/0012-9615(2006)076[0521:TTSDAW]2.0.CO;2 . Google Scholar CrossRef Search ADS Oliver , C.D. and Larson , B.C. 1996 Forest Stand Dynamics . John Wiley and Sons . Pompa-García , M. and Venegas-González , A. 2016 Temporal variation of wood density and carbon in two elevational sites of Pinus cooperi in relation to climate response in northern Mexico . PLoS One 11 , e0156782 . doi:10.1371/journal.pone.0156782 . Google Scholar CrossRef Search ADS PubMed Pretzsch , H. 2014 Canopy space filling and tree crown morphology in mixed-species stands compared with monocultures . For. Ecol. Manag. 327 , 251 – 264 . doi:10.1016/j.foreco.2014.04.027 . Google Scholar CrossRef Search ADS Pretzsch , H. and Biber , P. 2005 A re-evaluation of Reineke’ s rule and stand density index . For. Sci. 51 , 304 – 320 . Pretzsch , H. and Biber , P. 2016 Tree species mixing can increase maximum stand density . Can. J. For. Res. 46 , 1179 – 1193 . doi:10.1139/cjfr-2015-0413 . Google Scholar CrossRef Search ADS Pretzsch , H. , del Río , M. , Ammer , C. , Avdagic , A. , Barbeito , I. , Bielak , K. , et al. 2015 Growth and yield of mixed versus pure stands of Scots pine (Pinus sylvestris L.) and European beech (Fagus sylvatica L.) analysed along a productivity gradient through Europe . Eur. J. For. Res. 134 , 927 – 947 . doi:10.1007/s10342-015-0900-4 . Google Scholar CrossRef Search ADS Pretzsch , H. and Mette , T. 2008 Linking stand-level self-thinning allometry to the tree-level leaf biomass allometry . Trees 22 , 611 – 622 . doi:10.1007/s00468-008-0231-x . Google Scholar CrossRef Search ADS Pretzsch , H. and Schütze , G. 2015 Effect of tree species mixing on the size structure, density, and yield of forest stands . Eur. J. For. Res. 135 , 1 – 22 . doi:10.1007/s10342-015-0913-z . Google Scholar CrossRef Search ADS R-Core Team , 2016 . R: a language and environment for statistical computing. Reineke , L.H. 1933 Perfecting a stand-density-index for even-age forests . J. Agric. Res. 46 , 627 – 638 . Reyes-Hernandez , V. , Comeau , P.G. and Bokalo , M. 2013 Static and dynamic maximum size–density relationships for mixed trembling aspen and white spruce stands in western Canada . For. Ecol. Manag. 289 , 300 – 311 . doi:10.1016/j.foreco.2012.09.042 . Google Scholar CrossRef Search ADS Valbuena , P. , Peso , C. and Bravo , F. 2008 Stand density management diagrams for two mediterranean pine species in Eastern Spain . Invest. Agrar. Sist. Rec. 17 , 97 – 104 . Google Scholar CrossRef Search ADS Valladares , F. and Niinemets , U. 2008 Shade tolerance, a key plant feature of complex nature and consequences . Ann. Rev. Ecol. Evol. Syst. 39 , 237 – 257 . Google Scholar CrossRef Search ADS Vincent , G. and Harja , D. 2008 Exploring ecological significance of tree crown plasticity through three-dimensional modelling . Ann. Bot. 101 , 1221 – 1231 . doi:10.1093/aob/mcm189 . Google Scholar CrossRef Search ADS PubMed Vospernik , S. and Sterba , H. 2014 Do competition-density rule and self-thinning rule agree? Ann. For. Sci. 72 , 379 – 390 . doi:10.1007/s13595-014-0433-x . Google Scholar CrossRef Search ADS Weiskittel , A. , Gould , P. and Temesgen , H. 2009 Sources of variation in the self-thinning boundary line for three species with varying levels of shade tolerance . For. Sci. 55 , 84 – 93 . Woodall , C.W. , Miles , P.D. and Vissage , J.S. 2005 Determining maximum stand density index in mixed species stands for strategic-scale stocking assessments . For. Ecol. Manag. 216 , 367 – 377 . doi:10.1016/j.foreco.2005.05.050 . Google Scholar CrossRef Search ADS Wright , I.J. , Ackerly , D.D. , Bongers , F. , Harms , K.E. , Ibarra-Manriquez , G. , Martinez-Ramos , M. , et al. 2007 Relationships among ecologically important dimensions of plant trait variation in seven neotropical forests . Ann. Bot. 99 , 1003 – 1015 . doi:10.1093/aob/mcl066 . Google Scholar CrossRef Search ADS PubMed Zeide , B. 1987 Analysis of the 3/2 power law of self-thinning . For. Sci. 33 , 517 – 537 . Zeide , B. 2004 Optimal stand density: a solution . Can. J. For. Res. 34 , 846 – 854 . doi:10.1139/X03-258 . Google Scholar CrossRef Search ADS Zeller , L. , Ammer , C. , Annighöfer , P. , Biber , P. , Marshall , J. , Schütze , G. , et al. 2017 Tree ring wood density of Scots pine and European beech lower in mixed-species stands compared with monocultures . For. Ecol. Manag. 400 , 363 – 374 . doi:10.1016/j.foreco.2017.06.018 . Google Scholar CrossRef Search ADS Zhang , L. , Bi , H. , Gove , J.H. and Heath , L.S. 2005 A comparison of alternative methods for estimating the self-thinning boundary line . Can. J. For. Res. 35 , 1507 – 1514 . doi:10.1139/x05-070 . Google Scholar CrossRef Search ADS Zhang , J. , Oliver , W.W. and Powers , R.F. 2013 Reevaluating the self-thinning boundary line for ponderosa pine (Pinus ponderosa) forests . Can. J. For. Res. 43 , 971 . doi:10.1139/cjfr-2013-0133 . © Institute of Chartered Foresters, 2018. All rights reserved. For Permissions, please e-mail: journals.permissions@oup.com. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)

Journal

Forestry: An International Journal Of Forest ResearchOxford University Press

Published: Mar 21, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off