INTRODUCTION

The structure and dynamics of mixed natural forests are more difficult to characterize than those from even-aged stands given their considerable variability (^{Landres et al. 1999}). Forest variability is often expressed as heterogeneity, defined as any form of environmental variation, physical or biotic, occurring in space and/or time (^{Kuuluvainen et al. 2002}). Multi-scale factors can also influence individual tree responses and they can include: light quantity and quality, interception of precipitation from foliage, branches and litter, and root uptake of water and nutrients, among others (^{Kuuluvainen et al. 2002}). This small-scale heterogeneity often influences tree-to-tree competition, which in turn generates increased stand variability.

One of the most commercially and ecologically important forests in the southern hemisphere is composed of species from the Nothofagaceae family, which are found in Australia, New Zealand, New Guinea, New Caledonia, Chile and Argentina (^{Veblen et al. 1996}). The South American *Nothofagus* species dominate temperate and sub-Antarctic forests of Chile and Argentina. In Chile, there are nine species, including three evergreen and six deciduous species (^{Donoso 1993}). ^{Donoso (1993)} identified 12 forest types for the Chilean vegetation, one of these forest types is called RO-RA-CO. Its name originates from their three main *Nothofagus* species found there: roble (*Nothofagus obliqua* (Mirb.) Oerst*.*), raulí (*N. alpina* (Poepp. Endl.) Oerst.) and coihue (*N. dombeyi* (Mirb.) Oerst.). RO-RA-CO is the fourth largest forest type in Chile covering 10.8 % of the nation’s total area (^{CONAF 2011}) presenting a timber production level of 45 % of the total native wood market (~34,000 m^{3}, ^{INFOR 2016}).

The RO-RA-CO forest type is a complex ecosystem and is found in a vegetation gradient spanning wide geographical ranges in both latitude and altitude (^{Donoso 1993}). This forest type is found in southern Chile from 36° S to 41° S latitude approximately from the Pacific coast to the Andes mountain range. This wide distribution leads to a diverse set of ecosystems and growing conditions (^{Donoso 1993}, ^{Donoso et al. 1993}). These forests largely belong to the temperate hyperoceanic bio-climate zone with a large extension from coastal zones, interior valleys and mountainous areas, with diverse climates (^{Luebert and Pliscoff 2006}).

Several studies have addressed the patterns of heterogeneity in these forests. For example, ^{Echeverría and Lara (2004)} determined that the variability in diameter growth in *N. alpina* and *N. obliqua* was associated mainly with climate and soil characteristics. Another study found that, at a continental scale, *N. obliqua* productivity depends mainly on climate, although at a regional scale, edaphic and topographic variables are more relevant (^{Thiers et al. 2008}). Similarly, ^{Esse et al. (2013}) found that the geographical distribution of *N. dombeyi* is explained, in decreasing order of importance, by annual average rainfall, thermal oscillation and soil drainage. Two critical factors that affect the variability of RO-RA-CO forests are: 1) the presence of pure or mixed forests, and 2) their regeneration strategy (^{Donoso 1993}). Related to these factors, first, extreme climate and then soils conditions tend to favor pure forests (^{Donoso 1993}).

Forest disturbance is an important factor in the ecology of these *Nothofagus* spp. forests, as it contributes to increase its variability. Large scale disturbances, such as fire or earthquakes, usually promote a vast regeneration event, which promotes even-aged strata or cohorts (^{Luebert and Pliscoff 2006}). However, ^{Donoso (1993}) points out differences in the life history strategies of the *Nothofagus* species, which seem to be adapted to large-scale disturbance sizes and intensities. For example, *N. dombeyi* and *N. obliqua* are more of an early dominant pioneer species preferring large scale disturbances; whereas *N. alpina*, due to its intermediate shade tolerance, favors more patch-scale dynamics (^{Donoso 1993}, ^{Chauchard and Sbrancia 2003}, ^{Donoso et al. 2018}).

Tree- and stand-level growth has always been an important element used to characterize most natural forests and plantations, particularly when defining growth zones, which are useful to specify relatively homogeneous geographical areas with similar productivity and silvicultural response. For example, for *N. alpina* and *N. obliqua*, ^{Donoso et al. (1993}) defined four growth zones, where the most productive ones corresponded to those with warm temperatures. ^{Esse et al. (2013}) proposed for *N. dombeyi* five growth zones, which were defined, in decreasing order, by climate, soil and topography. Furthermore, ^{Gezan and Moreno (1999}) defined four growth zones and 11 sub-zones based on the inventory and evaluation of native vegetation resources of Chile (^{CONAF et al. 1999}), soil maps (^{Schlatter et al. 1994}, 1995), Agroclimatic map of Chile (^{INIA 1989}) and the growth zones proposed by ^{Donoso et al. (1993)}. Several studies on volume equations, basal area growth, dominant height-site and taper equations have found that the ^{Gezan and Moreno (1999)} zoning is useful to improve model predictability (^{Moreno 2001}, ^{Gezan et al. 2009}, ^{Moreno et al. 2017}).

One approach for better understanding the variability of any forest, and particularly the RO-RA-CO forest type, is by using multivariate analysis; which is a group of statistical methods that helps to describe and/or explore complex data as well as permitting formal inferences about the population (^{Everitt and Hothorn 2011}). One classification of multivariate analysis considers supervised and unsupervised methods (^{James et al. 2013}). In supervised methods, the goal is to predict a desired response from a set of predictor variables; whereas for unsupervised methods, the goal is to discover relationships among the data from the complete set of variables. For example, with the unsupervised methods of ordination or dimensionality reduction it is possible to understand data variability, or by using cluster analysis unknown similarities among objects can be discovered (^{James et al. 2013}). Some multivariate methods of ordination are principal component analysis (PCA), principal coordinates analysis (PCoA, also called multidimensional scaling), non-metric multidimensional scaling (NMDS), and correspondence analysis (CA) (^{Everitt and Hothorn 2011}). These methods help to reduce the complexity of data, visualize multidimensional data and develop research hypotheses (^{Everitt and Hothorn 2011}). Cluster analysis is a method that allows finding similar subgroups of objects in a dataset, where some of its commonly used variants are agglomerative hierarchical, K-means clustering, and model-based clustering (^{Everitt and Hothorn 2011}).

The aim of this study is to understand the underlying dasometric variability of second-growth forests of *N. obliqua, N. alpina* and *N. dombeyi*, through multivariate analysis of ordination (non-metric multidimensional scaling, principal coordinates analysis and principal components analysis) and clustering (K-means) using a group of tree- and stand-level variables typically obtained from forest inventories. In particular, we used forest inventory data from 158 plots (1,108 trees with individual diameter growth data) covering the complete geographical distribution of the RORACO forest type. Specific objectives were: 1) to reduce the dimensionality of data through ordination techniques; 2) to identify those variables or factors that explain high amount of data variability; and 3) to generate clusters to assess the growth zones proposed by ^{Gezan and Moreno (1999}).

METHODS

*Data description*

A total of 1,108 *N. obliqua*, *N. alpina* and *N. dombeyi* trees with diameter growth information were available for this analysis. These originated from a sample of 158 plots from a study by Universidad Austral de Chile (^{Ortega and Gezan 1998}) consisting of two different plot networks. Data were collected under a population stratification based on a national forest inventory (^{CONAF et al. 1999}). Further details of this dataset can be found in ^{Gezan and Moreno (1999)}.

This study used the average annual increment in diameter at breast height based on the penultimate and antepenultimate growth years, which was obtained from increment cores (864 trees) and stem sections (244 trees). Only two years were used to describe individual tree growth in order to maintain a strong relation with the current stand conditions. The last growth year was discarded as plots were established at different times within the growing season. The dataset includes several tree- and stand-level variables (table 1). The tree-level attributes considered were: species, diameter at breast height (DBH, cm), total height (H, m), breast height age (A, year), basal area of larger trees (BAL, m^{2} ha^{-1}), basal area of *Nothofagus* spp. of larger trees (BALn, m^{2} ha^{-1}), sociological status (SS, defined according to vertical stratification with dominant (1), codominant (2), intermediate (3) or suppressed (4)), relative basal area of larger trees (BALr, with BALr = BAL/BA) and the average annual increment in DBH (AIDBH, mm year^{-1}). Stand-level variables considered were: growth zone (Zone), stand basal area (BA, m^{2} ha^{-1}), density or number of trees per hectare (N, trees ha^{-1}), quadratic mean diameter (QD, cm), dominant height (Hd, m), dominant breast height stand age (Ad, years), dominant species (DS), site index (SI, m), basal area of *Nothofagus* spp. (BAN, m^{2} ha^{-1}), stand density index (SDI, trees ha^{-1}) and relative spacing (RS). Zones correspond to those proposed by ^{Gezan and Moreno (1999}), which are mainly based on climate, soil and potential site productivity (see appendix A for further details). Hd and Ad were calculated using the 100 largest trees per hectare in terms of DBH, and DS corresponded to the *Nothofagus* species with the largest proportion of BAN. SI was estimated using available dominant height-site models obtained on the same dataset for an index age of 20 years (see appendix B). BAN was calculated as the sum of basal area of *N. obliqua,*
*N. alpina,* and *N. dombeyi,* while SDI was calculated using: SDI = Nˣ(25.4/QD)^{β} (^{Avery and Burkhart 2002}), where β = -1.4112, as reported by ^{Gezan et al. (2007)}. Finally, RS was calculated according to the expression: RS = [(10,000/N)^{0.5}]/Hd (^{Avery and Burkhart 2002}).

Note: n, number of observations; SD, standard deviation; DBH, diameter at breast height (cm); H, total height (m); A, breast height age (years); BAL, basal area of larger trees (m^{2} ha^{-1}); BALn, basal area of larger trees for *Nothofagus* spp. (m^{2} ha^{-1}); SS, sociological status; BALr, relative BAL; AIDBH, annual increment in DBH (mm year^{-1}); BA, stand basal area (m^{2} ha^{-1}); N, number of trees (trees ha^{-1}); QD, quadratic diameter (cm); Hd, dominant height (m); Ad, dominant breast height stand age (years); SI, site index (m); BAN, basal area of *Nothofagus* spp. (m^{2} ha^{-1}); SDI, stand density index (trees ha^{-1}); RS, relative spacing.

*Multivariate methods*

To understand the underlying variability of the RORACO forest, three methods of ordination were evaluated: 1) non-metric multidimensional scaling (NMDS), 2) principal coordinate analysis (PCoA) and 3) principal component analysis (PCA). The goal of any ordination method is to represent main data trends through orthogonal axes (linearly independent, uncorrelated), thus reducing the original number of dimensions. The observed trends can then be interpreted graphically or further evaluated with the use of other multivariate techniques such as clustering (^{Borcard et al. 2011}).

NMDS is probably the most flexible ordination method that allows the use of continuous and categorical variables (^{Everitt and Hothorn 2011}), and unlike PCoA and PCA, it is not an eigenvalue-eigenvector-based method. NMDS, for the new data representation, keeps the order of the relationships among the initial set of objects but does not provide their exact multivariate distances (^{Borcard et al. 2011}). This method tolerates missing values, if there is enough information to identify the position of each object regarding others (^{Borcard et al. 2011}). For this study, the reduction of dimensionality was established using two axes (k = 2) with a modified Gower dissimilarity matrix (^{Anderson et al. 2006}). To compare variables at different scales, two NMDS analyses were performed: one with tree-level variables, and another with stand-level variables. Based on these analyses, ordination plots were used to represent the relationships among objects (trees or plots), where additional factors such as growth zones, tree species and dominant species were used in the identification of objects to help with the visualization of results. Here, object clouds were delimited using the extreme objects positions in relation to each factor level.

The PCoA, in contrast to NMDS, orders the data in axes according to their contribution to total variability based on eigenvalues, and allows for the use of any similarity/dissimilarity matrix, which represents the relationship among objects or variables. For the evaluation of PCoA in this study, both tree- and stand-level variables were analyzed jointly, where the number of coordinates were selected by examining a scree plot. Here, the same modified Gower dissimilarity matrix and graphical output generated for NMDS were also used.

The PCA is the main eigenvalue-eigenvector-based method that describes the variation of correlated quantitative variables, in terms of a new set of uncorrelated variables or axes. As with PCoA, these new axes or principal components (PC) are derived in decreasing order of importance according to their total variance explained (^{Everitt and Hothron 2011}). A PCA uses a similarity/dissimilarity matrix calculated based on the variances and covariances of the variables, or the correlations among them (^{Borcard et al. 2011}); thus, PCA uses the Euclidean distance to represent objects. Representation of objects (component scores) and original variables (factor loadings) can be done through a PCA Biplot. In this study, an exploratory PCA was used to identify relevant continuous variables. Here, two analyses were performed: one used both tree- and stand-level variables, and a second used only stand-level variables with the inclusion of the plot average AIDBH. For this multivariate method, to achieve normality assumptions, some variables were transformed using the inverse, logarithm, square root, or arc-sin functions. After transformations, variables were standardized. The number of PCA axes was determined by using (^{Quinn and Keough 2002}): 1) a latent root criterion by keeping components with eigenvalues > 1; 2) through a scree plot; and 3) by having a percent variance explaining over 90 %. Likewise, to facilitate the interpretation of the relationship among objects, distance-Biplot were used (^{Quinn and Keough 2002}).

K-means clustering is a method for dividing a dataset into K different and non-overlapping clusters. In this study, this multivariate method was implemented only with standardized stand-level variables. The determination of the number of clusters was done through scree plots of: 1) within-cluster sum of squares; and 2) by the silhouette method (^{Everitt and Hothorn 2011}). In addition, a Multiple Response Permutation Procedure (MRPP) was performed to test statistical differences among clusters (^{McCune et al. 2002}). The Pearson’s chi-squared test of association was used to establish any association between clusters and growth zones (as defined by ^{Gezan and Moreno 1999}). If any relationship was found, a post-hoc test for all pairwise comparisons with Bonferroni corrections of P-values was performed (^{MacDonald and Gardner 2000}).

All statistical analyses were performed using R v. 3.2.2 (^{R Core Team 2015}). For all analyses, a significance level of 5 % was used. The “vegan” version 2.4-1 R package was used for dissimilarities matrices, NMDS, PCoA, and MRPP (^{Oksanen 2011}). “Cluster” version 2.0-4 was used for the silhouette plot (^{Maechler 2013}).

RESULTS

Many of the analyzed tree- and stand-level variables presented, as expected, a high level of correlation between them, such as BA-BAN, Hd-RS, BA-SDI, BAL-BALn, Hd-QD, A-Ad, RS-QD and DBH-H, with values > 0.83; thus, indicating a high level of multicollinearity.

The multivariate method of NMDS, using categorical and continuous data for tree-level variables, converged with a stress goodness-of-fit value of 0.09, corresponding to a good representation of the original distance matrix. Stress values below 0.10 are considered reasonable, and values below 0.05 are often an indication of a good fit. For this group of variables, NMDS ordination with the two selected axes did not identify different groups of objects, only a single group is represented, which is surrounded by isolated data points (figure 1). For this representation, there were no clear patterns for the different growth zones or species, *i.e.* growth zones represented almost exactly the same tree-level cloud, and *N*. *obliqua* had only minor differences when compared to *N. alpina* and *N. dombeyi*. This indicates that these factors do not explain much of the variability. Nevertheless, the axis 1 is related to tree-to-tree competition variables, specifically BAL, BALr, BALn and SS; and the axis 2 has variables associated with tree size and growth, such as A, H, DBH and AIDBH.

Another NMDS analysis was performed only with the stand-level variables (figure 2), resulting in a stress value of 0.05, which also represents a good description of the original distance matrix. As with the tree-level variables, the clouds for growth zones were overlapping; however, the dominant species showed a small difference between *N.*
*obliqua* and *N*. *dombeyi*, while *N*. *alpina* was practically contained within the other two dominant species clouds. Here, the first axis can be described as “stand density” because of the strong effect of N. The second axis is related to variables associated with stand development, such as BA, QD, BAN, SDI, Hd and Ad. Interestingly, none of the productivity variables (SI, Zone) indicated relevant importance in explaining data variability.

The ordination method of PCoA explained ~90 % of the cumulative variance using two dimensions (PCoA axes) when both tree- and stand-level variables were evaluated together (data not shown). As expected, we observed that trees belonging to a plot were very close (overlap) in these new coordinates (figure 3); hence, it appears that the variability associated with the stand-level variables dominates over the one with individual tree variables. In addition, the results from PCoA confirms the findings from NMDS, with overlapping growth zones and larger similarities between *N. alpina* and *N. dombeyi*. Furthermore, N is a relevant variable explaining ~78.5 % of the variability (axis 1, figure 3). Moreover, the axis 2 is described by stand development stage variables (as with NMDS for the stand-level subset, figure 2), and several tree-level variables (A, H, BAL, BALn, specie (Sp), AIDBH) however with less relevance. Finally, productivity variables, as with NMDS analysis, showed a low level of importance in explaining variability.

The third ordination method PCA was evaluated for both combined groups of variables (tree- and stand-level), to assist in identification of variables that better explained the variability of the sample population. To achieve normality, the variables DBH, BA, BALn and AIDBH were transformed using a square root, and BALr was transformed with the arcsine function. The results of PCA indicated that four or five principal components (PC) explain most of the variability (88 % and 92 %, respectively) in these RORACO forests. Factor loadings for these five PC are presented in table 2. PC1 (with 37 % of the variance explained) can be described as a “stand development stage” axis; PC2 (24 %) a “tree-to-tree competition” axis; PC3 (16 %) a “stocking and stand competition” axis; PC4 (11 %) a “site productivity” axis; and PC5 (4 %) a “tree DBH growth” axis.

Note: DBH, diameter at breast height (cm); H, total height (m); A, breast height age (years); BAL, basal area of larger trees (m^{2} ha^{-1}); BALn, basal area of larger trees for *Nothofagus* spp. (m^{2} ha^{-1}); SS, sociological status; BALr, relative BAL; AIDBH, annual increment in DBH (mm year^{-1}); BA, stand basal area (m^{2} ha^{-1}); N, number of trees (trees ha^{-1}); QD, quadratic diameter (cm); Hd, dominant height (m); Ad, dominant breast height stand age (years); SI, site index (m); BAN, basal area of *Nothofagus* spp. (m^{2} ha^{-1}); SDI, stand density index (trees ha^{-1}); RS, relative spacing.

PC1 is described mostly by the stand-level variables BA, QD, Hd, Ad, BAN, SDI and RS. The other stand-level variables N and SI were relevant only on the third and fourth components, respectively. In relation to tree-level variables, we observed a high contribution in PC1 for those related to tree size (DBH, H, and A), while tree-to-tree competition variables (BAL, BALn and BALr) were important in PC2 with almost identical factor loadings. AIDBH was also strongly associated with this component but with opposite direction (*i.e.*, negative factor loading). A Biplot for PC1 and PC2 explained 61 % of the variability, and is presented in figure 4A. Recall that PC1 is mainly representing stand development stages, and PC2 tree-to-tree competition; hence, the component scores associated with each observation show a wide range of tree-to-tree competition at later stand development stages (left side of figure 4A). For early stand development stages, this tree-to-tree competition is more uniform (right side of figure 4A). This Biplot also shows the negative association between AIDBH and tree competition variables; where, as expected, high individual competition means low diameter growth. Surprisingly, productivity (represented by SI) had a negligible effect in explaining the variability of this population. From this Biplot, groups of correlated variables can be easily identified: Hd-Ad-QD, BA-BAN-RS, BAL-BALn-BALr-AIDBH and A-N. In addition, no clear patterns of the distribution of the trees could be attributed to growth zone and species.

A second PCA was done with only stand-level variables and average plot AIDBH. Factor loadings of the four components are shown in table 3. PCs, as expected, were almost identical to the ones described earlier when both groups of variables were used, with PC1 (43 % variance explained) as a “stand development stage” axis; PC2 (24 %) a “stocking and stand competition” axis; PC3 (18 %) a “site productivity and growth” axis; and PC4 (8 %) as “stand competition and growth” axis. The two first PCs explained 67 % of the total variance and they are represented by a Biplot in figure 4B. Here, plot observations (component scores) present a higher variability in stand competition at later stand development stages (right side of figure 4B), and *vice versa*.

Note: DBH, diameter at breast height (cm); H, total height (m); A, breast height age (years); BAL, basal area of larger trees (m^{2} ha^{-1}); BALn, basal area of larger trees for *Nothofagus* spp. (m^{2} ha^{-1}); SS, sociological status; BALr, relative BAL; AIDBH, annual increment in DBH (mm year^{-1}); BA, stand basal area (m^{2} ha^{-1}); N, number of trees (trees ha^{-1}); QD, quadratic diameter (cm); Hd, dominant height (m); Ad, dominant breast height stand age (years); SI, site index (m); BAN, basal area of *Nothofagus* spp. (m^{2} ha^{-1}); SDI, stand density index (trees ha^{-1}); RS, relative spacing.

The implementation of K-means clustering with stand-level variables identified two clusters using within-cluster sum of squares and silhouette methods. The first cluster was associated with high values of N and RS, while the second with high values of BA, BAN, Ad, Hd and QD (figure 4B and table 4). Hence, the main differences among plots were associated with stand development stages, where, first, there are young forests with high tree density, and second, older forests with basal area concentrated in a few individuals. An MRPP test statistic of ‘chance corrected within-group agreement’ (As) of 0.16 was obtained for these two clusters indicating strong statistical differences between them. In addition, the observed delta statistic was 3.51 (*P* < 0.001) versus an expected delta value of 4.17, also indicating strong significant differences (*P* < 0.001) in terms of stand-level variables. Association between clusters and growth zones was analyzed using a Pearson’s chi-squared statistic, that resulted in a value of 8.69 (*P* = 0.034). Post-hoc, pairwise comparisons only established a small association between clusters and the growth zones 1 and 2 (*P* = 0.049). These two clusters are presented in a map together with the growth zones, where no clear spatial patterns were noted between them (figure 5).

Note: AIDBH, annual increment in DBH (mm year^{-1}); BA, stand basal area (m^{2} ha^{-1}); N, number of trees (trees ha^{-1}); QD, quadratic diameter (cm); Hd, dominant height (m); Ad, dominant breast height stand age (years); SI, site index (m); BAN, basal area of *Nothofagus* spp. (m^{2} ha^{-1}); SDI, stand density index (trees ha^{-1}); RS, relative spacing.

DISCUSSION

This study analyzed the variability of second-growth RO-RA-CO forests of *N. obliqua*, *N. alpina*, and *N. dombeyi* forests, using multivariate analysis of ordination, clustering and multi-scale typically measured in forest inventories. Overall, we found high levels of variability in tree- and stand-level variables, with wide ranges across growth zones in terms of averages and standard deviation (table 1). This variability is likely due to several factors including biophysical attributes, auto-ecology and species-specific characteristics, as well as disturbance type and level, among other factors (^{Donoso 1993}, ^{Chauchard and Sbrancia 2003}, ^{Echeverria and Lara 2004}, ^{Luebert and Pliscoff 2006}).

In Chile, the Andean mountain vegetation is marked by altitudinal gradients that are affected by temperature and rainfall, where higher altitudes are characterized by lower temperatures and increased rainfall (^{Luebert and Pliscoff 2006}). Additionally, the latitudinal gradient strongly affects vegetative growth, southern forests have a shorter period (^{Donoso 1993}). To simplify the effect of biophysical variables on a forest resource and its silvicultural management responses, it is useful to define relatively homogeneous geographical growth zones. The zoning used in this study (^{Gezan and Moreno 1999}) was defined by combining different sources of information, such as climate, soil, vegetation distribution and site productivity. In this study, we evaluated this zoning using more advanced statistical techniques such as multivariate analysis of ordination and clustering.

Results show that the unsupervised multivariate methods identified important variables at the tree- and stand-level for the RO-RA-CO forest type in Chile. Both PCAs performed, using only continuous variables, exhibited and explained a high percentage of variance, ranging from 37 % to 43 %. At the stand-level the most important variables were N, BA, BAN, QD, Ad and Hd. In addition, two clear clusters were found, which represent two conditions of stand development stages: 1) young forests with a high tree density; and 2) adult forests with high basal area concentrated in few large trees. Other authors have reported differences based on stand developmental stages. For example, ^{Donoso (1993}) found, for *Nothofagus* spp., changes in diameter distribution from unimodal to bimodal as the stand ages, reflecting additional heterogeneity that is present in these forests. This heterogeneity increases with the interaction with biophysical variables such as climate, altitude and latitude.

The other multivariate methods, NMDS and PCoA, using categorical and/or continuous data, evidenced a small differentiation among species, although not for growth zones. ^{Lusk and Ortega (2003}) found no differences in productivity for comparisons based on the same zoning. The results of NMDS and PCoA at the stand level showed that, in the multivariate space, *N*. *alpina* was intermediate between the other two species, being closer to *N. dombeyi* than to *N. obliqua* in relation to dasometric variables from our dataset. Other authors also have found important differences among these species. For example, ^{Chauchard and Sbrancia (2003}) studied shade tolerance, and their results indicated that *N. alpina* is the most tolerant followed by *N. dombeyi*, and then *N. obliqua*. In terms of tree-to-tree competition, the same authors reported that at older ages, *N. dombeyi* stands present higher density than that presented by *N. obliqua*, which is an indication of the differential levels of stocking between these species, the latter also was indicated by ^{Lusk and Ortega (2003)}. ^{Donoso (1993}) documents that *N*. *obliqua* has a strong natural thinning at 40-50 years of age due to its low ability to tolerate tree-to-tree competition. Similar results regarding the importance of stocking were found in this study, where the NMDS and PCoA analysis showed that N was the most relevant variable to explain variability. In addition, our results show that tree-to-tree competition, as described by PC2, explained 24 % of the total variability for both tree- and stand-level variables, and this component was mainly described by the variables BAL, BALr and BALn. Interestingly, these variables were also important for NMDS, although not for PCoA analyses. Another relevant result associated with tree-to-tree competition is the strong relationship that exists between BALr and SS, where the first variable can be used to completely replace the more subjective field observed sociological status (see figure 1)

Several studies have identified a wide range of factors that influence the variability found in the RO-RA-CO forest type, as were exposed early. Biophysical factors were not directly included in this study, nonetheless growth zones essentially function as proxies for such factors. We admit that this is a limitation when trying to explain forest variability only with tree- and stand-level variables. However, we note that this study focused on the understanding of the behavior of these dasometric variables as an initial step to assist with silviculture planning and on the future development of growth and yield model for RO-RA-CO forest in Chile.

In conclusion, this study contributed to a better understanding of the underlying variability for an important forest type in Chile using multivariate statistical analyses based on tree- and stand-level variables typically obtained from traditional forest inventories. Our findings show, for the datasets evaluated, the limited role that site productivity and the current growth zones have in explaining the variability of this sampled population. Conversely, we note that stand development stages (N, BA, BAN, QD, Hd and Ad), tree-to-tree competition (BAL, BALr, BALn, and SS) and size-tree (DBH, H and A) attributes are the critical variables needed to explain this variability. Hence, the identification of these factors, through the use of multivariate tools, can be critical to understand current silvicultural responses and to plan future management regimes of these forests, and they are essential elements to consider in future endeavors of modelling these complex mixed forests.