SciELO - Scientific Electronic Library Online

 
vol.39 número3Parámetros estructurales de la vegetación arbóreo-arbustiva del bioma Caatinga sometida a sistemas silviculturales en la región semiárida de BrasilDiversidad genética del cedro rojo ( Cedrela odorata ) en el estado de Tabasco, México índice de autoresíndice de materiabúsqueda de artículos
Home Pagelista alfabética de revistas  

Servicios Personalizados

Revista

Articulo

Indicadores

Links relacionados

  • En proceso de indezaciónCitado por Google
  • No hay articulos similaresSimilares en SciELO
  • En proceso de indezaciónSimilares en Google

Compartir


Bosque (Valdivia)

versión On-line ISSN 0717-9200

Bosque (Valdivia) vol.39 no.3 Valdivia  2018

http://dx.doi.org/10.4067/S0717-92002018000300397 

Artículos

Exploring stand and tree variability in mixed Nothofagus second-growth forests through multivariate analyses

Explorando la variabilidad de rodal y árbol en bosques de renovales mixtos de Nothofagus spp. utilizando análisis multivariados

Paulo C Moreno a  

Salvador A Gezan b  

Sebastian Palmas b  

Francisco J Escobedo c  

Wendell P Cropper Jr. b  

a Centro de Investigación en Ecosistemas de la Patagonia, Coyhaique, Chile

b University of Florida, School Forest Resources & Conservation, 363 Newins-Ziegler Hall, PO Box 110410, Gainesville, FL, 32611, USA, phone: +1 (352) 846.0133

c Universidad del Rosario, Faculty of natural Sciences & Mathematics, Bogota, Colombia

SUMMARY:

Second-growth forests of Nothofagus obliqua (roble), N. alpina (raulí) and N. dombeyi (coihue), known locally as RO-RA-CO forest type, are among the most important natural mixed forest types of Chile. Several studies have identified a wide range of factors that could influence both stand and tree variability found in these forests. To better characterize potential tree- and stand-level factors that are associated with RO-RA-CO variability, and that are available in typical forest inventories, several unsupervised multivariate statistical methods were evaluated: 1) non-metric multidimensional scaling (NMDS); 2) principal coordinates analysis (PCoA); and 3) principal component analysis (PCA). The data used in this study originated from a sample of 158 plots consisting of two plot networks that covered the full geographic area of the RO-RA-CO forest type in Chile. We found that site productivity and growth zones did not explain the differences within the sampled population. However, stand development stages, tree-to-tree competition, and tree-size attributes were critical variables with a high percentage of variance explained using PCA, ranging from 61 % to 67 %. In addition, for the PCoA analysis, the variable stand density is important, with ~78 % variance explained.

Key words: RO-RA-CO forest type; PCA; NMDS; PCoA; K-means cluster

RESUMEN:

Los bosques de renovales de Nothofagus obliqua (roble), N. alpina (raulí), y N. dombeyi (coihue), conocidos como el tipo forestal RORACO, es uno de los más importantes bosques nativos mixtos de Chile. Varios estudios han identificado un amplio rango de factores que influyen en la variabilidad a nivel de rodal y árbol encontrada en estos bosques. Este estudio, contribuye a una mejor caracterización potencial de los factores de árbol y rodal que están asociados con la variabilidad de los bosques de RORACO, y que son obtenidas desde inventarios forestales tradicionales. Varios métodos estadísticos no supervisados de análisis multivariado fueron evaluados: 1) análisis de escalado multidimensional no métrico (NMDS), 2) análisis de coordenadas principales (PCoA) y 3) análisis de componentes principales (PCA). La base de datos de este estudio está sustentada en una muestra de 158 parcelas formada a partir de dos redes de parcelas para toda el área geográfica del tipo forestal RORACO en Chile. Nuestros resultados indican que la productividad de sitio y zonas de crecimiento no explican considerablemente las diferencias encontradas en esta población muestreada. Sin embargo, el estado de desarrollo del rodal, la competencia individual árbol a árbol, y los atributos de tamaño de árbol, son variables críticas, con un alto porcentaje de la varianza explicada usando PCA, alcanzando un rango de entre 61 % a 67 %. Además, para el análisis de PCoA, la variable de densidad del rodal es importante con un ~78 % de la varianza explicada.

Palabras clave: tipo forestal RORACO; métodos multivariados; componentes principales; clúster de K-promedios

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 m3, 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, m2 ha-1), basal area of Nothofagus spp. of larger trees (BALn, m2 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, m2 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, m2 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).

Table 1 Summary statistics of this study dataset for Nothofagus obliqua, N. dombeyi and N. alpina tree- and stand-level variables separated by growth zones according to Gezan and Moreno (1999). Resumen de estadísticos de la base de datos de estudio para variables de árbol y rodal de Nothofagus obliqua, N. alpina, y N. dombeyi separado por zona de crecimiento definidas de acuerdo a Gezán y Moreno (1999).  

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 (m2 ha-1); BALn, basal area of larger trees for Nothofagus spp. (m2 ha-1); SS, sociological status; BALr, relative BAL; AIDBH, annual increment in DBH (mm year-1); BA, stand basal area (m2 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. (m2 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.

Figure 1 Non-Metric Multidimensional Analysis (NMDS) using the Nothofagus obliqua, N. alpina, and N. dombeyi tree-level variables, identifying growth zones (A) and species (B). Análisis de escalado multidimensional no métrico (NMDS) usando la base de árbol de Nothofagus obliqua, N. alpina, y N. dombeyi, identificando zonas de crecimiento (A) y especies (B). 

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.

Figure 2 Non-Metric Multidimensional Analysis (NMDS) using the Nothofagus obliqua, N. alpina, and N. dombeyi stand-level variables, identifying growth zones (A) and species (B). Análisis de Escalado Multidimensional No Métrico (NMDS) usando la base de rodal de Nothofagus obliqua, N. alpina, y N. dombeyi, identificando zonas de crecimiento (A) y especie dominante (B). 

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.

Figure 3 Principal coordinates analysis (PCoA) using the Nothofagus obliqua, N. alpina, and N. dombeyi tree- and stand-level variables together, identifying growth zones, ZONE (A) and species, SPd (B). Análisis de Coordenadas Principales (PCoA) usando las bases de árbol y rodal de Nothofagus obliqua, N. alpina, y N. dombeyi combinadas, identificando zonas de crecimiento, ZONE (A) y especies, SPd (B). 

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.

Table 2 Factor loading for the first five principal components (PC1 to PC5) obtained from the Nothofagus obliqua, N. alpina and N. dombeyi tree- and stand-level variables together, with the transformation used indicated. Values in parenthesis correspond to the proportion of the variance explained by each component. Factor loadings smaller than 0.1 in absolute value are not shown, and those larger than 0.3, in absolute value, are in bold.Cinco primeros componentes principales (PC1 a PC5) obtenidos de las dos bases de datos de árbol y rodal combinadas de Nothofagus obliqua, N. alpina, y N. dombeyi, junto a la transformación usada en el análisis. Valores en paréntesis corresponden a la proporción de la varianza explicada por cada componente. Valores menores a 0,1, en valor absoluto, no se muestran, y aquellos mayores a 0,3, en valor absoluto, están en negrita.  

Note: DBH, diameter at breast height (cm); H, total height (m); A, breast height age (years); BAL, basal area of larger trees (m2 ha-1); BALn, basal area of larger trees for Nothofagus spp. (m2 ha-1); SS, sociological status; BALr, relative BAL; AIDBH, annual increment in DBH (mm year-1); BA, stand basal area (m2 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. (m2 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.

Figure 4 First two PCA axes with scores and factor loadings obtained using tree- and stand-level variables (A) and stand-level variables together with average plot AIDBH (B) of Nothofagus obliqua, N. alpina, and N. dombeyi. The two K-means clusters are also identified in (B).Primeros dos ejes de componentes principales (PCA) con objetos y valores propios obtenidos usando: la base de árbol y rodal combinadas (A), y la base de rodal más el promedio por parcela de AIDBH (B) para Nothofagus obliqua, N. alpina, y N. dombeyi. Los dos conglomerados obtenidos por el método K-promedios están identificados en (B). 

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.

Table 3 Factor loading for the first four principal component (PC1 to PC4) obtained from the Nothofagus obliqua, N. alpina and N. dombeyi stand-level variables together with average plot AIDBH, with the used transformation indicated. Values in parenthesis correspond to the proportion of the variance explained by each component. Factor loadings smaller than 0.1 in absolute value are not shown, and those larger than 0.3, in absolute value, are in bold. Cuatro primeros componentes principales (PC1 a PC4) obtenidos de las dos bases de datos de rodal de Nothofagus obliqua, N. alpina, y N. dombeyi más el promedio por parcela de AIDBH combinadas, junto a la transformación usada en el análisis. Valores en paréntesis corresponden a la proporción de la varianza explicada por cada componente. Valores menores a 0,1, en valor absoluto, no se muestran, y aquellos mayores a 0,3, en valor absoluto, están en negrita.  

Note: DBH, diameter at breast height (cm); H, total height (m); A, breast height age (years); BAL, basal area of larger trees (m2 ha-1); BALn, basal area of larger trees for Nothofagus spp. (m2 ha-1); SS, sociological status; BALr, relative BAL; AIDBH, annual increment in DBH (mm year-1); BA, stand basal area (m2 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. (m2 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).

Table 4 Average values of clusters obtained from Nothofagus obliqua, N. alpina and N. dombeyi stand-level variables together with average plot AIDBH for the K-means clustering method. Valores promedios de conglomerados obtenido de la base de rodal de Nothofagus obliqua, N. alpina, y N. dombeyi más el promedio por parcela de AIDBH usando el método de agrupación de K-promedios. 

Note: AIDBH, annual increment in DBH (mm year-1); BA, stand basal area (m2 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. (m2 ha-1); SDI, stand density index (trees ha-1); RS, relative spacing.

Figure 5 Geographical location of Nothofagus obliqua, N. alpine and N. dombeyi sampled plots identified according to the cluster from the K-means clustering method. Growth zones defined by Gezan and Moreno (1999) are also included.Localización geográfica de las parcelas de Nothofagus obliqua, N. alpina, y N. dombeyi muestreadas identificadas de acuerdo a los conglomerados del método de agrupación K-promedios. Las zonas de crecimiento definidas por Gezan y Moreno (1999) también están identificadas. 

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.

ACKNOWLEDGEMENTS

This research was supported by University of Florida, through the School of Forest Resources & Conservation. The authors acknowledge to Dr. Alicia Ortega and Dr. Guillermo Trincado from Universidad Austral de Chile for facilitating the two Nothofagus plot datasets.

REFERENCES

Anderson MJ, KE Ellingsen, BH McArdle. 2006. Multivariate dispersion as a measure of beta diversity. Ecology Letters 9:683-693. [ Links ]

Avery TE, HE Burkhart. 2002. Forest measurements, 5th edn. New York, USA. McGraw-Hill. 456 p. [ Links ]

Borcard D, Gillet F, Legendre P. 2011. Numerical ecology with R. New York, USA. Springer. 306 p. [ Links ]

Chauchard L, R Sbrancia. 2003. Diameter growth models for Nothofagus obliqua. Bosque 24(3):3-16. [ Links ]

CONAF (Corporación Nacional Forestal, CL), CONAMA (Comisión Nacional del Medio Ambiente, CL), BIRF (Banco Internacional de Reconstrucción y Fomento, USA), Universidad Austral de Chile, Pontificia Universidad Católica de Chile, Universidad Católica de Temuco. 1999. Proyecto catastro y evaluación de los recursos vegetacionales nativos de Chile. Santiago, Chile. Ministerio de Agricultura. 88 p. [ Links ]

CONAF (Corporación Nacional Forestal, CL). 2011. Catastro de los recursos vegetacionales nativos de Chile. Monitoreo de cambios y actualizaciones. Periodo 1997 - 2011. Santiago, Chile. Ministerio de Agricultura. 25 p. [ Links ]

Donoso C. 1993. Bosques templados de Chile y Argentina. Variación, estructura y dinámica. Santiago, Chile. Editorial Universitaria. 483 p. [ Links ]

Donoso PJ, C Donoso, V Sandoval. 1993. Proposición de zonas de crecimiento de renovales de roble (Nothofagus obliqua) y raulí (Nothofagus alpina) en su rango de distribución natural. Bosque 14(2):37-55. [ Links ]

Donoso PJ, A Promis, D Soto. 2018. Silvicultura en bosques nativos. Experiencias en silvicultura y restauración en Chile, Argentina y el oeste de Estados Unidos. Valdivia, Chile. Chile Inititiatve, OSU College of Forestry, America. 260 p. [ Links ]

Echeverria C, A Lara. 2004. Growth patterns of secondary Nothofagus obliqua-N. alpina forests in southern Chile. Forest Ecology and Management 195:29-43. [ Links ]

Esse C, PJ Donoso, V Gerding, F Encina-Montoya. 2013. Determination of homogeneous edaphoclimatic zones for the secondary forests of Nothofagus dombeyi in central-southern Chile. Ciencia e Investigación Agraria 40(2):351-360. [ Links ]

Everitt B, T Hothorn. 2011. An introduction to applied multivariate analysis with R. New York, USA. Springer. 273 p. [ Links ]

Gezan SA, PC Moreno. 1999. Establecimiento y medición de una red de parcelas permanentes para renovales de roble (Nothofagus obliqua), raulí (Nothofagus alpina) y coigüe (Nothofagus dombeyi). Bosque Nativo 22:9-11. [ Links ]

Gezan SA, PC Moreno, A Ortega. 2009. Modelos fustales para renovales de roble, raulí y coigüe en Chile. Bosque 30(2):61-69. [ Links ]

Gezan SA, A Ortega, E Andenmatten. 2007. Diagramas de manejo de densidad para renovales de roble, raulí y coigüe en Chile. Bosque 28(2):97-105. [ Links ]

INIA (Instituto de Investigación Agropecuaria, CL). 1989. Mapa agroclimático de Chile. Proyecto agrometereología. Santiago, Chile. Ministerio de Agricultura. [ Links ]

INFOR (Instituto Forestal, CL). 2016. Bosque nativo 12. Ministerio de Agricultura. Accessed 13 March 2017. Available in Available in http://wef.infor.cl/publicaciones/bn/2016/12/BN201612.pdf . [ Links ]

James G, D Witten, T Hastie, R Tibshirani. 2013. An introduction to statistical learning with applications in R. New York, USA. Springer. 426 p. [ Links ]

Kuuluvainen T. 2002. Natural variability of forests as a reference for restoring and managing biological diversity in boreal Fennoscandia. Silva Fennica 36(1):97-125. [ Links ]

Landres PB, P Morgan, FJ Swanson. 1999. Overview of the use of natural variability concepts in managing ecological systems. Ecological Applications 9:1179-1188. [ Links ]

Luebert F, P Pliscoff. 2006. Sinopsis bioclimática y vegetacional de Chile. Santiago, Chile. Editorial Universitaria. 316 p. [ Links ]

Lusk CH, A Ortega. 2003. Vertical structure and basal area development in second-growth Nothofagus stands in Chile. Journal of Applied Ecology 50:639-645. [ Links ]

MacDonald PL, RC Gardner. 2000. Type I error rate comparisons of post hoc procedures for I J chi-square tables. Educational and Psychological Measurement 60(5):735-754. [ Links ]

Maechler M. 2013. Cluster analysis extended Rousseeuw et al. R CRAN. [ Links ]

McCune B, JB Grace, DL Urban. 2002. Analysis of ecological communities. MjM Software Design, Gleneden Beach. Oregon, USA. [ Links ]

Moreno PC. 2001. Proposición preliminar de curvas de índice de sitio para renovales de raulí. Tesis Ingeniero Forestal. Santiago, Chile. Facultad de Ciencias Forestales, Universidad de Chile. 49 p. [ Links ]

Moreno PC, S Palmas, FJ Escobedo, WP Cropper, SA Gezan. 2017. Individual-tree diameter growth models for mixed Nothofagus second-growth forests in southern Chile. Forests 8(12):506-525. [ Links ]

Oksanen J. 2011. Multivariate analysis of ecological communities in R: vegan tutorial. R package version 1.17-7 [ Links ]

Ortega A, SA Gezan. 1998. Cuantificación de crecimiento y proyección de calidad en Nothofagus. Bosque 19(1):123-126. [ Links ]

Quinn GP, MJ Keough. 2002. Experimental design and data analysis for biologists. Cambridge, United Kingdom. Cambridge University Press. 537 p. [ Links ]

R Core Team. 2015. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/Links ]

Schlatter J, V Gerding, J Adriazola. 1994. Sistema de ordenamiento de la tierra. Herramienta para la planificación forestal, aplicada a las regiones VII, VIII y IX. Valdivia, Chile. Serie técnica, Facultad Ciencias Forestales, Universidad Austral de Chile. 114 p. [ Links ]

Schlatter J, V Gerding, H Huber. 1995. Sistema de ordenamiento de la tierra. Herramienta para la planificación forestal, aplicada para la X región. Valdivia, Chile. Serie técnica, Facultad Ciencias Forestales, Universidad Austral de Chile. 110 p. [ Links ]

Thiers O, V Gerding, E Hildebrand. 2008. Renovales de Nothofagus obliqua en centro y sur de Chile: Factores de sitio relevantes para su productividad. In Bava J, Picco OA, Pildaín MB, López P, Orellana I eds. Libro de actas de Eco Reuniones: segunda reunión sobre los Nothofagus en la Patagonia. Esquel, Argentina. CIEFAP. p. 255-260. [ Links ]

Veblen TT, RS Hill, J Read. 1996. The ecology and biogeography of Nothofagus forests. New Haven, USA. Yale University Press. 414 p. [ Links ]

1Gezan and Moreno (2000). Curvas de Sitio – Altura Dominante para Renovales de Roble, Raulí y Coigue. Documento Interno. Proyecto FONDEF D97I1065. Universidad Austral de Chile.

APPENDIX A: RO-RA-CO growth zones proposed by Gezan and Moreno (1999)

This zoning was built in three stages: 1) boundary delimitation of RO-RA-CO forests, 2) analysis of biophysical variables, and 3) growth zones definition. The delimitation in stage one was based on the cartography of inventory and evaluation of native vegetation resources of Chile (CONAF et al. 1999), selecting those second-growth forest stands dominated by Nothofagus obliqua, N. alpina and N. dombeyi. This delimitation contained different stands corresponding to others land uses; thus, incorporating not only the actual land use but potential areas for these species.

The second stage incorporated cartography information related to soil, climate, geomorphology, and vegetation. The studies used were the following:

  1. Soil maps (Schlatter et al. 1994, 1995). These studies proposed a zoning based on Köppen climate classification (macroscale) and soil classes (microscale).

  2. Agro-climatic map of Chile (INIA 1989). Corresponding to a cartography of homogeneous climatic zones based on soil variables.

  3. Growth zones (Donoso et al. 1993). This study proposed a zoning for second-growth forests, based on different field studies and expert opinion.

The third and final stage combined all above information and delimited potential growth zones and subzones following these main criteria:

  1. Separate the mild temperate climate (Cfa) and warm-summer Mediterranean climate (Csb) with the oceanic climate (Cf) using the Köppen climate classification.

  2. Altitudinal classification using the soil maps creating three categories: 0-400, 400-800, and > 800 m a.s.l.

  3. Latitudinal classification using geomorphology changes, such as rivers and climatic data.

  4. The growth zones proposed by Donoso et al. (1993) were used as a framework for N. obliqua and N. alpina considering the similar delimitations for zones and subzones.

  5. The agro-climatic map was used to delimit subzones on those growth zones representing large geographical areas or contrasting soil types.

The classification results showed four zones and eleven subzones (see figure 5) with the following broad characteristics:

Zone 1. It occupies the western slopes of the Coast Range in Bio-Bio region and the central valley between Biobío region and Valdivia province. This zone includes the Mediterranean climate with oceanic influence, but with dry summer most years. It includes all forests of N. alpina and N. dombeyi of the Nahuelbuta Range, in other areas N. obliqua is predominant. It is divided into three subzones.

Zone 2. It corresponds entirely to the Andean Range lower slopes between Bio-Bio region and Valdivia province. It is characterized by a high rainfall but with low relative air humidity in summer and the presence of frosts. Nothofagus obliqua and N. dombeyi stands are common in this zone. It includes two subzones.

Zone 3. It includes the westerns slope of the Andean Range between Araucania region to Llanquihue province, excluding the higher areas of this mountain range. High rainfall is common but with a short dry period (< two months). It is dominated by N. obliqua stands but N. dombeyi appears in wet areas. Four growth subzones conform this zone.

Zone 4. It corresponds to high Andean Range, and their western boundary coincides with the N. alpina presence. The highest rainfall and lowest temperatures are typical in this zone. The three species are present, including the majority of N. alpina stands of the Andean Range between Biobío to Los Lagos regions. It includes two subzones.

APPENDIX B: Dominant height-site model

The dominant height-site models were developed by Gezan and Moreno (2000)1, using the data that originated from the project reported by Ortega and Gezan (1998). This model has the following general equation:

where, Hd = dominant height (m); Ad = dominant age (years); SI = site index to 20 years (m); and a, β 0 , β 1 are the estimated parameters (table B-1). This was fitted with the same database, target species, and growth zones of this study.

Table B-1 Parameter estimates of the dominant height-site model by specie and growth zone. Parámetros estimados clasificados por especie y zona de crecimiento para los modelos de indicie de sitio - altura dominante. 

Received: August 23, 2017; Accepted: July 10, 2018

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License