Regional Variability of the Romanian Main Tree Species Growth Using National Forest Inventory Increment Cores

In many countries, National Forest Inventory (NFI) data is used to assess the variability of forest growth across the country. The identification of areas with similar growths provides the foundation for development of regional models. The objective of the present study is to identify areas with similar diameter and basal area growth using increment cores acquired by the NFI for the three main Romanian species: Norway spruce (Picea abies L. Karst), European beech (Fagus sylvatica L.), and Sessile oak (Quercus petraea (Matt.) Liebl.). We used 6536 increment cores with ages less than 100 years, a total of 427,635 rings. The country was divided in 21 non-overlapping ecoregions based on geomorphology, soil, geology and spatial contiguousness. Mixed models and multivariate analyses were used to assess the differences in annual dimeter at breast height and basal area growth among ecoregions. Irrespective of the species, the mixed models analysis revealed significant differences in growth between the ecoregions. However, some ecoregions were similar in terms of growth and could be aggregated. Multivariate analysis reinforced the difference between ecoregions and showed no temporal grouping for spruce and beech. Sessile oak growth was separated not only by ecoregions, but also by time, with some ecoregions being more prone to draught. Our study showed that countries of median size, such as Romania, could exhibit significant spatial differences in forest growth. Therefore, countrywide growth models incorporate too much variability to be considered operationally feasible. Furthermore, it is difficult to justify the current growth and yield models as a legal binding planning tool.


Introduction
The development of national forest growth and yield models served the society needs of the last century, but fails to address the current challenges posed by the rise in population, expected standard of life, and constant changes of the environment [1]. The perceived endless computational power available today recommends the distillation of the national models, such that they are tailored to specific questions or areas. Such detailed models address not only the society's expectations from the forest, but also serve as the basis of sustainable forest management. Arguably, the next refining modeling level is the identification of areas for which distinct models can be developed. The natural identification of such areas are the ecoregions. A plethora of definitions are available for an ecoregion, Figure 1. The PSC of the Romanian NFI overlaid on the general location of the forest, in green. The color of the PSC reflects the geomorphology: red for plains, blue for hills, and black for mountains. The spatial distribution of the PSCs does not reflect the actual density [12].
The Romanian NFI collected a large number of increment cores, which represent valuable information for the development, refinement and spatial delineation of regional growth models. Analyzing NFI growth data in combination with regional spatial units, such as the ecoregions, represents the next logical step in the development of improved growth and yield models. Therefore, the objective of this study is to identify the existence of spatial variation in tree growth for the three main forest species in Romania: Norway spruce (Picea abies (L) H. Karst), European beech (Fagus sylvatica L), and Sessile oak (Quercus petraea (Matt.) Liebl). Furthermore, if the presence of spatial variation is confirmed, different patterns of annual tree growth are to be examined based on the characteristics of the ecoregions. Even though other studies have used increment cores from inventory data to analyze forest growth [13][14][15], the use of NFI increment cores data on a large area for model development represents, to our knowledge, a novelty in the literature.

Romanian Ecoregions
One of the main drivers that could be responsible for tree growth variation across the country, if differences in growth exist, are likely the ecoregions. In this study, an ecoregion is defined as a contiguous area with similar geology, soils and geomorphology. The information describing the drivers of an ecoregion was supplied by the Romanian Academy. Based on the four drivers, we divided the country into 21 ecoregions (Figure 2). The ecoregions vary in size from 161,904 ha, the The PSC of the Romanian NFI overlaid on the general location of the forest, in green. The color of the PSC reflects the geomorphology: red for plains, blue for hills, and black for mountains. The spatial distribution of the PSCs does not reflect the actual density [12].
The Romanian NFI collected a large number of increment cores, which represent valuable information for the development, refinement and spatial delineation of regional growth models. Analyzing NFI growth data in combination with regional spatial units, such as the ecoregions, represents the next logical step in the development of improved growth and yield models. Therefore, the objective of this study is to identify the existence of spatial variation in tree growth for the three main forest species in Romania: Norway spruce (Picea abies (L.) H. Karst), European beech (Fagus sylvatica L.), and Sessile oak (Quercus petraea (Matt.) Liebl). Furthermore, if the presence of spatial variation is confirmed, different patterns of annual tree growth are to be examined based on the characteristics of the ecoregions. Even though other studies have used increment cores from inventory data to analyze forest growth [13][14][15], the use of NFI increment cores data on a large area for model development represents, to our knowledge, a novelty in the literature.

Romanian Ecoregions
One of the main drivers that could be responsible for tree growth variation across the country, if differences in growth exist, are likely the ecoregions. In this study, an ecoregion is defined as a contiguous area with similar geology, soils and geomorphology. The information describing the drivers of an ecoregion was supplied by the Romanian Academy. Based on the four drivers, we divided the country into 21 ecoregions (Figure 2

Increment Cores
The species considered in this study, namely Norway spruce, European beech and Sessile oak, not only cover approximately 60% of the Romanian forests, but are also the main species in terms of economic, ecological, and social importance. To identify the presence of spatial variation in the growth of the tree species, we have used the increment cores of the Romanian NFI. The cores are extracted at each 5-year measurement at breast height (1.3 m) parallel with the contour line. The number of trees cored from each SP depends on the number of species: when only one species was present, then 3-4 cores were extracted, otherwise 2-3 cores/species were collected. For each cored tree, only one increment core was extracted and the cored trees are selected randomly, though are conditioned by specific constraints. The trees must belong to the dominant and co-dominant canopy class [17] and should not exhibit exterior damage. The cores are further mounted on a solid support, sanded, scanned and the width of each ring is measured using a graphical procedure. The main information extracted from the increment cores is the width of each ring and their succession. To ensure validity of the comparisons, we have included only the increment cores that were extracted from even-aged pure species stands. In addition, to ensure compatibility among various ecoregions, only stands with no thinning or no active management were considered. A more detailed description of the data can be found in Marin et al. [16].
The data used in this study consists of 6536 increment cores from 1655 PSC the last one being collected in 2011. The total number of rings is 427,635, out of which 136,904 rings are for Norway spruce, 241,240 rings for European beech, and 49,491 rings for Sessile oak. To focus the study on ages that are operational relevant, only rings with maximum age of 99 years were included. The number

Increment Cores
The species considered in this study, namely Norway spruce, European beech and Sessile oak, not only cover approximately 60% of the Romanian forests, but are also the main species in terms of economic, ecological, and social importance. To identify the presence of spatial variation in the growth of the tree species, we have used the increment cores of the Romanian NFI. The cores are extracted at each 5-year measurement at breast height (1.3 m) parallel with the contour line. The number of trees cored from each SP depends on the number of species: when only one species was present, then 3-4 cores were extracted, otherwise 2-3 cores/species were collected. For each cored tree, only one increment core was extracted and the cored trees are selected randomly, though are conditioned by specific constraints. The trees must belong to the dominant and co-dominant canopy class [17] and should not exhibit exterior damage. The cores are further mounted on a solid support, sanded, scanned and the width of each ring is measured using a graphical procedure. The main information extracted from the increment cores is the width of each ring and their succession. To ensure validity of the comparisons, we have included only the increment cores that were extracted from even-aged pure species stands. In addition, to ensure compatibility among various ecoregions, only stands with no thinning or no active management were considered. A more detailed description of the data can be found in Marin et al. [16].
The data used in this study consists of 6536 increment cores from 1655 PSC the last one being collected in 2011. The total number of rings is 427,635, out of which 136,904 rings are for Norway spruce, 241,240 rings for European beech, and 49,491 rings for Sessile oak. To focus the study on ages that are operational relevant, only rings with maximum age of 99 years were included. The number of rings (cores) for each ecoregion varies according to the species (Table 2), from 272 rings (4 cores) (i.e., European beech in the Maramures Plateau (323)) to 49,576 (986 cores) (i.e., Norway spruce in the Eastern Carpathians (131)). Each ring was described with five attributes: ecoregion, species, year, age, and width. The year is the calendar year of each ring within the increment cores, age represents the age of each ring, and ring width (RW) is the ring width in micrometers. Diameter at breast height (DBH) and basal area (BA) increments are derived from ring width, and are considered more appropriate to characterize tree growth than RW. To avoid values for DBH growth determined by special environmental conditions, we have used only cores from trees located inside the natural distribution range. Therefore, for Norway spruce, we have only included increment cores from the mountainous regions in the analysis, for European beech only those from the hill, plateau and mountainous regions, and for Sessile oak only those from the hill and plateau regions.

Analysis
To assess the differences in annual growth among ecoregions, we have used mixed model and multivariate procedures. Each type of procedure addresses a different facet of the relationship between ecoregions and growth: the mixed model reflects the annual growth within a neighborhood of one year, whereas the multivariate considers the annual growth within the context of all years. In the absence of meteorological data, we implicitly considered the change in climate by constraining the analysis to the last 50 calendar years (i.e., since 1960). We have chosen as threshold the year 1960, as multiple studies pointed a shift in climate occurring at the end of the six decade of the twentieth century, not only in Europe [18][19][20] but also in the world [21][22][23]. Climate is an additional major driver of tree growth, besides site, but in this study we focus only on site, which forces the usage of a narrow time interval for assessment of spatial variability.

Mixed Models Analysis
To assess the impact of the ecoregion on the radial growth at breast height we have used a linear mixed modelling approach covariates age and DBH for DBH increment, or age and BA for BA  (1) and (2)). The ecoregions were treated as fixed effects whereas the age was considered as random effect, as the age or each ring depends on the year.
We introduced age in the model to account for difference in radial growth with age [24][25][26]. We selected DBH and BA as covariate to account for the impact of the size of the tree on the growth.
Given the temporal dependency of the measurements and the nested structure of the data, we implemented a repeated measurements framework within the liner mixed model. Ecoregion and the covariate were considered as a fixed effect, whereas age was considered random. Among the possible covariance structures, we have considered the compound symmetry, variance components, and autoregressive of order 1 [27,28]. To identify the appropriate structure we have used the Akaike's Information Criterion (AIC) [29]. To formally determine the differences among ecoregions, if any, we have used the estimated marginal means, as recommended by Gotway and Stroup [30]. To detect the existence of residual autocorrelation, we used the Durbin-Watson test [31]. All analyses were performed in SAS 9.4 [32].

Multivariate Analysis
To simultaneously examine the spatial growth across all ecoregions we have used multivariate analysis having as variables width of the annual rings. Therefore, the multivariate analysis should have for each individual tree, and consequently each ecoregion, a series of at least 10 variables, depending on its age. However, trees grow differently through time, which means that rings with the same age must also be grouped by year. The requirement of having rings with the same age for each ecoregion and year since 1960 was not fulfilled. Therefore, instead of using individual trees, we have used groups of trees that have the same age and year. Even with this consolidation of data, old rings will be rare, meaning that only few trees would be considered in analysis. Consequently, the ring age should be reduced to a value that ensures at least one tree per calendar year.
Among the multivariate methods available for the analysis of the ring widths, we have used principal component analysis (PCA), hierarchical cluster analysis, and canonical discrimination analysis. The PCA represents an orthogonal transformation of the variable matrix that is used for variable reduction and data exploration [32]. We performed PCA to investigate if there is any evidence of grouping the growth by ecoregion or if this grouping is determined by a subset of ring widths. Hierarchical cluster analysis has a similar objective as PCA because it hierarchically clusters the observations in the dataset, but using an unsupervised method as the number of clusters is not established apriori. We have created clusters using the Ward method [33]. The cubic clustering criterion and the scree plot analysis of eigenvalues were employed in determining the number of clusters [34,35]. We complemented the hierarchical cluster analysis with canonical discrimination analysis, which grouped the ecoregions starting with an a priori number of groups [36]. We selected the number of groups of ecoregions as the number of clusters identified by the cluster analysis. The difference among the ecoregions was tested with Wilks' Lambda, Pillai's Trace, the Hotelling-Lawley Trace and Roy's Greatest Root [33]. We executed all the analyses in SAS 9.4 [32].

Mixed Models Analysis
The tree growth charts revealed an obvious impact of the ecoregion on the diameter increment and basal area increment ( Figure 3). Depending on the species and the variable measuring growth, there are ecoregions that clearly have different growth increments than the rest, such as the Eastern Carpathians (131) for Norway spruce, the Maramures plateau (323) for European beech, or the Caras  (221), and the Transylvania Plateau (324) for the DBHi of Sessile oak. The inference on ecoregion impact on growth based on visual information was supported by the summary statistics for each species and attribute, as the increments varies by more than 100% for DBHi and BAi for European beech and Sessile oak ( Table 3). The descriptive statistics paint a compelling argument, considering that the variance does not necessarily vary among ecoregions ( Table 3). The summary statistics and the charts suggest different increments for some ecoregions, which support the formal analysis of the differences among ecoregions with mixed models. 8 inference on ecoregion impact on growth based on visual information was supported by the summary statistics for each species and attribute, as the increments varies by more than 100% for DBHi and BAi for European beech and Sessile oak ( Table 3). The descriptive statistics paint a compelling argument, considering that the variance does not necessarily vary among ecoregions ( Table 3). The summary statistics and the charts suggest different increments for some ecoregions, which support the formal analysis of the differences among ecoregions with mixed models.   While visual inspection and the summary statistics supported the hypothesis of the study that there are differences in increments among ecoregions, its formal assessment could reveal the existence of a multitude of different growth increments. Therefore, the inclusion of covariates in the analysis seemed justified, as the visual inspection of the relationships DBHi-DBH and BAi-BA revealed not only differences among some ecoregions, but also a trend, particularly for BAi (Figure 4). While visual inspection and the summary statistics supported the hypothesis of the study that there are differences in increments among ecoregions, its formal assessment could reveal the existence of a multitude of different growth increments. Therefore, the inclusion of covariates in the analysis seemed justified, as the visual inspection of the relationships DBHi-DBH and BAi-BA revealed not only differences among some ecoregions, but also a trend, particularly for BAi ( Figure  4). Irrespective of the species, the AIC identified the first order autoregressive structure as the most suitable covariance structure. For Norway spruce, the ecoregions were found to have a significant effect on annual DBH growth and basal area growth. The highest estimate for average DBH increment was found in the Banat Mountains (331) ecoregion, and the lowest in the West Southern Carpathians (233) ecoregion. The estimated marginal means showed that most of the ecoregions are significantly different from one another in terms of growth (Table 4), except for eight pairs. In total, we found evidence that the DBH growth of Norway spruce significantly differs in six distinct ecoregions: the Irrespective of the species, the AIC identified the first order autoregressive structure as the most suitable covariance structure. For Norway spruce, the ecoregions were found to have a significant effect on annual DBH growth and basal area growth. The highest estimate for average DBH increment was found in the Banat Mountains (331) ecoregion, and the lowest in the West Southern Carpathians (233) ecoregion. The estimated marginal means showed that most of the ecoregions are significantly different from one another in terms of growth (Table 4), except for eight pairs. In total, we found evidence that the DBH growth of Norway spruce significantly differs in six distinct ecoregions: the Eastern Carpathians (131), Buzau-Vrancea mountains (231), East Southern Carpathians (232), West Southern Carpathians (233), Banat Mountains (331), and Volcanic Ridge (333). The seventh region analyzed, the Western Carpathians (332), was found to be similar with four other regions. There is no evidence that the results are artifacts [37], as the main assumption needed to be checked in repeated measurements analysis, namely that residuals are white noise, was fulfilled (i.e., the Durbin-Watson test indicated no autocorrelation). The growth of European beech also varies significantly by ecoregion, the highest average DBH increment being found in the Maramures Plateau (323), and the lowest in the Caras Hills (321) ( Table 3). However, not all regions were significantly different from one another, mirroring the Norway spruce. In the case of European beech, we noticed groups of ecoregions that have both similar DBH and BA increments (Table 5) (Table 5). Similarly to Norway spruce, there was no evidence that the results are based on methods violating analytical assumptions, as the Durbin-Watson test indicated no autocorrelation of the residuals.
Sessile oak mirrored the findings for Norway spruce and European beech, the ecoregions having a significant effect on tree growth ( Table 6). The highest average growth was in the Moldavian Hills (122) ecoregion and the lowest in the Caras Hills (321) ( Table 3). The regions exhibited mostly different growth from one another, with the exception of three pairs: the Buzau-Vrancea piedmonts (221) and Caras Hills (321); Caras Hills (321) and Cris Hills (322); Caras Hills (321) and Transylvania Plateau (324) ( Table 6). There are also other four other similar pairs, but only from a DBH increment perspective. The Durbin-Watson test also found no autocorrelation of the residuals. It can be concluded that there are seven ecoregions spanning the Romanian areal of Sessile oak that significantly differ in tree growth: the Moldavian Plateau (121), Moldavian Hills (122), Buzau-Vrancea piedmonts (221), Getic Plateau (222), Cris Hills (322), Maramures plateau (323) and Transylvania Plateau (324). Caras Hills (321) was found similar in DBH or BA increment with four other ecoregions.

Multivariate Analysis
As expected, a value less than 50 ensured at least one tree per calendar year and ecoregion. Depending on the species, the number of annual rings that was used in all multivariate analyses is either 20 for sessile oak or 30 for Norway Spruce and European beech. Irrespective of the attribute measuring the growth, DBHi or BAi, the three multivariate analyses supported the finding of mixed model investigation for all three species, as a significant separation among regions was revealed. The PCA, for which the first two principal components explained at most 80% of the growth variation, which occurred for BAi of sessile oak, (i.e., 79% for Norway Spruce and 73% for European beech), did not reveal an evident grouping according to the ecoregion ( Figure 5). Therefore, globally, there are other aspects of tree growth which are more important in delineating growth than ecoregions.

Multivariate Analysis
As expected, a value less than 50 ensured at least one tree per calendar year and ecoregion. Depending on the species, the number of annual rings that was used in all multivariate analyses is either 20 for sessile oak or 30 for Norway Spruce and European beech. Irrespective of the attribute measuring the growth, DBHi or BAi, the three multivariate analyses supported the finding of mixed model investigation for all three species, as a significant separation among regions was revealed. The PCA, for which the first two principal components explained at most 80% of the growth variation, which occurred for BAi of sessile oak, (i.e., 79% for Norway Spruce and 73% for European beech), did not reveal an evident grouping according to the ecoregion ( Figure 5). Therefore, globally, there are other aspects of tree growth which are more important in delineating growth than ecoregions. According to the cubic clustering criterion, the hierarchical cluster analysis found three groups of ecoregions for all three species (Figure 6), regardless the attribute measuring growth. Mirroring PCA findings, the groups were not as clearly delineated as mixed models analysis, but the results show distinct growing patterns among the three areas as pointed by the minimum difference among the groups ( Table 7). The clear separation of the ecoregion is the results of the amount of information covered by the first two or three eigenvectors, which are responsible for the groups, that is at least 78%. Furthermore, the difference among the groups was more than 15% (the case of Norway spruce), but almost 50% for European beech and Sessile oak (i.e., >45%). The hierarchical clustering revealed that Norway spruce growing in the ecoregions Eastern Carpathians (131), Western Carpathians (332) and the volcanic ridge (333) are similar and distinct from the rest. The results were less conclusive for the European beech and sessile oak, as a clear pattern did not emerge from the tree hierarchy. According to the cubic clustering criterion, the hierarchical cluster analysis found three groups of ecoregions for all three species (Figure 6), regardless the attribute measuring growth. Mirroring PCA findings, the groups were not as clearly delineated as mixed models analysis, but the results show distinct growing patterns among the three areas as pointed by the minimum difference among the groups ( Table 7). The clear separation of the ecoregion is the results of the amount of information covered by the first two or three eigenvectors, which are responsible for the groups, that is at least 78%. Furthermore, the difference among the groups was more than 15% (the case of Norway spruce), but almost 50% for European beech and Sessile oak (i.e., >45%). The hierarchical clustering revealed that Norway spruce growing in the ecoregions Eastern Carpathians (131), Western Carpathians (332) and the volcanic ridge (333) are similar and distinct from the rest. The results were less conclusive for the European beech and sessile oak, as a clear pattern did not emerge from the tree hierarchy.  Sessile oak 3 0.25 71.5 87.5 Among all the multivariate methods, the canonical discrimination analysis was the only one that showed an obvious difference in growth between ecoregions (Figure 7), regardless of the attribute (i.e., DBHi or BAi). All four multivariate tests indicated that for all three species, the class means vectors are not equal. It was clear that Norway spruce from the Western Carpathians (131) grow differently than the West Southern Carpathians (233), and both from the rest of the ecoregions ( Figure  7a). Whereas PCA and hierarchical cluster analysis failed to identify ecoregions with distinct growth from European beech, the canonical discriminant analysis succeeded, as the Moldavian plateau (121) and Moldavian hills (122) exhibit different growing patterns from the rest and between them ( Figure  7b). We also found a clear delineation of ecoregions by growth for sessile oak, with the Buzau-Vrancea Piedmonts (221) showing different DBHi and BAi than the Transylvania plateau (324), which are both distinct from the other ecoregions. It should be noticed that fewer ecoregions can be considered by the canonical discriminant analysis (Figure 7c), as the ecoregions with less trees, consequently less factorial combinations, are not included (e.g., in the case of sessile oak, which has only four ecoregions, compared with eight in the mixed model analysis).   Among all the multivariate methods, the canonical discrimination analysis was the only one that showed an obvious difference in growth between ecoregions (Figure 7), regardless of the attribute (i.e., DBHi or BAi). All four multivariate tests indicated that for all three species, the class means vectors are not equal. It was clear that Norway spruce from the Western Carpathians (131) grow differently than the West Southern Carpathians (233), and both from the rest of the ecoregions (Figure 7a). Whereas PCA and hierarchical cluster analysis failed to identify ecoregions with distinct growth from European beech, the canonical discriminant analysis succeeded, as the Moldavian plateau (121) and Moldavian hills (122) exhibit different growing patterns from the rest and between them (Figure 7b). We also found a clear delineation of ecoregions by growth for sessile oak, with the Buzau-Vrancea Piedmonts (221) showing different DBHi and BAi than the Transylvania plateau (324), which are both distinct from the other ecoregions. It should be noticed that fewer ecoregions can be considered by the canonical discriminant analysis (Figure 7c), as the ecoregions with less trees, consequently less factorial combinations, are not included (e.g., in the case of sessile oak, which has only four ecoregions, compared with eight in the mixed model analysis).  Sessile oak 3 0.25 71.5 87.5 Among all the multivariate methods, the canonical discrimination analysis was the only one that showed an obvious difference in growth between ecoregions (Figure 7), regardless of the attribute (i.e., DBHi or BAi). All four multivariate tests indicated that for all three species, the class means vectors are not equal. It was clear that Norway spruce from the Western Carpathians (131) grow differently than the West Southern Carpathians (233), and both from the rest of the ecoregions ( Figure  7a). Whereas PCA and hierarchical cluster analysis failed to identify ecoregions with distinct growth from European beech, the canonical discriminant analysis succeeded, as the Moldavian plateau (121) and Moldavian hills (122) exhibit different growing patterns from the rest and between them ( Figure  7b). We also found a clear delineation of ecoregions by growth for sessile oak, with the Buzau-Vrancea Piedmonts (221) showing different DBHi and BAi than the Transylvania plateau (324), which are both distinct from the other ecoregions. It should be noticed that fewer ecoregions can be considered by the canonical discriminant analysis (Figure 7c), as the ecoregions with less trees, consequently less factorial combinations, are not included (e.g., in the case of sessile oak, which has only four ecoregions, compared with eight in the mixed model analysis).

Discussion
Tree growth is influenced by a wide range of factors across a multitude of gradients [38]. Knowledge about this variation can lead to improved and localized forest growth models, with direct impact on forest management, forest operations and planning [39][40][41]. Our analysis was able to identify spatial/regional variations in the annual growth of the main tree species even for a country close to the median size (ranked 81 by IndexMundi), such as Romania.
The mixed models analysis showed that growth significantly differs by ecoregion. For all three species, several ecoregions grouped, indicating that regionalization of various forest related attributes or models is required. In the case of Norway spruce, there was similarity in DBH and BA increment between ecoregions Western Southern Carpathians (233) and Banat Mountains (331). One explanation for the analogous growth can be geomorphology, as all regions were mountainous, have similar soils and aspect (south facing) and are spatially adjacent. Spatial proximity can represent an important factor when comparing ecoregions because it can encompass numerous factors that influence growth. Among the most cited of those factors is climate, with effects at multiple levels in tree physiology [42][43][44], but others such as genetics, disturbances or localized stress events can also have an impact on tree development [45]. The DBH growth in the Banat Mountains (331) was also not different from Western Carpathians (332) and Volcanic ridge (333) ecoregions. This correspondence could be based on geomorphological factors, as soils, geology and location are different or it can also be attributed to other confounding factors. DBH growth also did not differ between Western Southern Carpathians (233) and Western Carpathians (332) and Volcanic ridge (333). As the ecoregions Banat Mountains (331) and Volcanic ridge (333) show similar DBH increment and mostly BA increment, these two ecoregions are good candidates for being joined from the perspective of Norway spruce growth. The highest average growth (Banat Mountains (331)) and the lowest (West Southern Carpathians (233)) occurred in a mountainous region with a similar soil composition, but a different altitude. Banat Mountains (331) are having a lower altitude with stepped leveling appearance and remnant limestone geology, whereas West Southern Carpathians (233) present higher altitudes, with leveled relief and limestone geology. The altitude can represent a proxy for other variables influencing growth, such as temperature, precipitation or for areas prone to stress and disturbances.
In the case of European beech, similar growth was found between pairs of plateau regions (ecoregions Buzau-Vrancea piedmonts (221), Getic plateau (222), Caras Hills (321), Transylvanian plateau (324)). As beech has its native range at these altitudes, it is expected that growth in these ecoregions is higher and also closer in terms of increments. However, the fact that these ecoregions are not completely equivalent, illustrates that the spatial variability of growth is present even within same geomorphology and that other factors can be involved in determining tree growth. These can be soil or geology, but also climatic variables or the presence of localized stress events. The average growth in the Maramures Plateau (323) is different from all the other ecoregions, with the highest basal area increment estimate. One possible explanation for this difference could be attributed to the combined effects of growing conditions and other factors such as regional climate or natural distribution of the genetic material. The growth similarity between Eastern Carpathians (131) and some plateau regions (Moldavian Plateau (121), Caras Hills (321), Transylvania Plateau (324)) indicates again that there are multiple factors influencing growth, and that geomorphology is not the only defining factor.
For Sessile oak, only a few ecoregions showed growth similarities. The Caras Hills (321) ecoregion exhibited similar growth with four other ecoregions where Sessile oak is present (Buzau-Vrancea piedmonts (221), Getic plateau (222), Cris Hills (322) and Transylvania Plateau (324)), all being located at the foothills of the mountains. This indicates that this foothill ecoregion represents a location of average growth when compared with other Sessile oak ecoregions, and could be combined with the others. Highest average growth occurred in a hilly ecoregion (Moldavian Hills (122)) and a much lower one in the neighboring one, the Moldavian Plateau (121). The unambiguous difference between these two regions underlines the spatial differentiation in tree growth for Sessile oak. As in Norway spruce, the altitude represents a potential differentiation factor, with the plateau ecoregion being around 200-500 m high and the hills around 500-900 m high.
Two models for tree growth were developed in this study, and ecoregions could be differentiated by both models. Some authors [46,47] suggest that both DBH and BA increment are appropriate in growth modelling, while Russell [48] suggests that DBH increment might decrease the model error more than BA. The results of this study showed that, for all three species and especially for European beech and Sessile oak, both the DBH increment and BA increment models converged towards same results in a high number of cases. This is a clear indication that ecoregions have an effect on tree growth, irrespective of the type of variable used to measure increment. There were instances where only DBH increment indicated a growth differentiation between ecoregions, as this variable might be more suitable for increment modelling [48].
The multivariate analysis examined simultaneously all the DBHi or BAi, looking at the growth pattern along the life of each individual tree. The PCA and cluster analysis suggested that DBHi cannot evidently group ecoregions and cannot be grouped by certain periods of time. This means that within the same species, variation of growth in time did not represent a grouping factor. The canonical discrimination analysis indicated that there is growth difference by ecoregion for both Norway spruce and European beech. However, for Sessile oak, the canonical variables have discriminatory power for classifying growth into ecoregions, irrespective of DBHi or BAi. As the canonical variables are a linear combination of yearly growth, DBHi in certain years could determine the inclusion in an ecoregion. Sessile oak has its native range in warm areas which are more prone to draught and having water as limiting factor. Although further analysis in combination with climate data is required, it is possible that certain ecoregions were more influenced by changes in climatic conditions, which were instrumental in this growth differentiation of ecoregions.
One main finding of the mixed models analysis was that geomorphology plays a major role in the spatial variation of growth. Ecoregions were found to have similar growth increments, even when soil, geology or location differ. The geomorphology was found to be relevant for all three species. One explanation for this finding is that geomorphology changes the specific soil characteristics [49] and indirectly influences growth. It is also possible that geomorphology has driven the historical human impact on forest cover and growth (i.e., higher impact in flatter areas and lower in mountains ones), and therefore influencing in this alternate manner the vegetation development [50]. Another explanation can be that climate conditions present in that geomorphology type have a strong impact on tree growth [44,51,52], which can be a stronger driver than soil [53]. Non-climatic variables such as geomorphology can represent a proxy for different climatic characteristics [54,55], as spatially adjacent ecoregions with different altitudes exhibited very different growth increments. However, even when considering only non-climate variables, spatial variability was evident.

Conclusions
The spatial variation of growth for tree species represents an essential aspect of forest management planning, modelling and reporting [56]. This study aimed at determining the existence of variation in DBH and BA growth over the area of a medium sized country located in a temperate region. A series of 21 ecoregions summarized the ecological conditions for the entire country.
Using mixed models and multivariate analysis on 6536 increment cores from the Romanian NFI, we were able to establish that tree growth is significantly different between various ecoregions for Norway spruce, European beech and Sessile oak. Numerous factors such as geology, soil and geomorphology are responsible for the difference in growth, with the latter being frequently present in the differentiation of ecoregions. An important finding of this study was that global, country wide, growth models incorporate far too much variability from an operational perspective. They can be used for country-wide resource estimations, but it is hard to justify their usage as a legal binding planning tool.
As the change in climate also influences the vegetation growth, and consequently, the type and amount of products and services supplied by each ecosystem, further dendrological studies that include climate and temporal analysis are needed for a better understanding of tree growth variability.
Author Contributions: G.M. organized data collection and framed the study, V.C.S. interpreted and wrote the manuscript, I.V.A. wrote portions of the manuscript, B.M.S. analyzed, interpreted, and wrote portions of the manuscript. All authors have read and agreed to the published version of the manuscript.