The Polish Provenances of European Larch Overperform the Expected Growth Dynamics Indicated by the Sigmoid Model

: This article attempts an unusual interpretation of the observations characteristic of experiments that compare different tree species’ provenances. The focus falls on larch ( Larix decidua subsp. polonica ). The data came from the experiment established in 1967 at the Siemianice forest experimental station, where 21 Polish larch provenances were compared. The study’s main objective was to compare the basic estimates of growth dynamics, the maximum growth rate and acceleration, and the point in time when these values were achieved. A four-parameter sigmoid growth function was used to model the average stand basal area increments and its ﬁrst and second derivatives to calculate the indicators of the growth dynamic. The models explained 98% of the observed variations resulting from the 21 inventories. Only one growth parameter showed a statistically signiﬁcant difference among the compared provenances. The G ó ra Chełmowa provenance achieved the highest value of the maximum growth acceleration, but it was statistically signiﬁcantly different only from the three underperformed provenances. However, when the average values for all the experimental plots ( n = 86) were compared with those of the deterministic model (the stand volume and yield tables), the maximum growth rate and acceleration values were higher for the former. We discuss the potential factors responsible for this overperformance and point out the potential risks that arise from growth and quality metrics only when deciding on the best-performing provenances. The sigmoid growth model employed in this study might be an excellent alternative for comparing the growth dynamics among different stands or replications in experimental studies. Considering only the early results, the sigmoid growth model proves its limitations, and the conclusions reached should be treated with caution.


Introduction
European larch (Larix decidua Mill.) is one of the four gymnosperms of high importance for forestry in Poland. The distribution range of larch was shaped after the last glacial period and currently encompasses the Alps, the Sudetes, and a scattered range in the Carpathians up to the Apuseni Mountains [1]. It is currently cultivated on a larger scale in the lowlands [2]. Altitudinally, it can grow up to the tree lines in the Alps. As a result of taxonomical debate, L. decidua is assigned to different subspecies. However, Pâques et al. [3] indicated only two subspecies: L. decidua subsp decidua Mill. and L. decidua subsp. polonica (Racib. ex Wóycicki) Domin. The latter species dominates in Poland, as stated by the cited authors.
The first attempts to introduce larch in the lowlands for cultivation purposes showed that the Alpine origins were characterized by slow growth, an early bud break, and thus became susceptible to late frosts when grown there [4]. The initial comparisons among the European provenances investigated by Schober [5] indicated that the provenances from the Sudetes, the Tatra Mountains, and those from the lower altitudes of the Alps were characterized by a relatively high growth rate. Furthermore, the susceptibility to Lachnellula willkommii (R. Hartig) Dennis, an agent of European larch canker, was much lower among Sudeten provenances than Alpine ones. The IUFRO initiative established two international experiments: the first in 1944 focused on 48 provenances (plus one for Larix kaempferi Lambert and one for Larix sibirica Ledebour), and the second one in 1958/69 focused on 63 provenances [4]. Since then, numerous national experiments have been established. Three aspects that are important in gathering knowledge in the field of forestry are time, space, and environmental variability. Hence, following Giertych's [6] viewpoint, comparing experiments conducted simultaneously and replicated across several other localities is of the greatest value. A study by Szymański and Wilczyński [7] is an excellent example of provenance experiment results compared among three different experimental sites. Suppose the circumstances mentioned above are met. In that case, another issue is to coordinate forest owners to manage and conduct observations at similar intervals and in the same manner to comply with the ceteris paribus rule.
Sigmoid functions play a crucial role in the modelling of tree growth dynamics [8]. It was broadly discussed which one performs best based on the goodness-of-fit and predictability. For instance, Colbert et al. [9] compared four different functions and indicated that the Chapman-Richards function was prominent. To date, two studies employed a similar approach to compare the differences in growth among provenances [10,11]. Both studies were carried out in the early stages of stand development: at 11 and 14 years old. The most significant advantage of modelling tree growth with sigmoid growth functions is the possibility of interpreting the function graph. The analysis of the first and second sigmoid model derivatives allows for the calculation of growth dynamics [8] rather than a static interpretation of a point in time, e.g., the volume of the stand when ready for harvesting. Here, we present the interim results of a single Polish experimental site-provenance trial. The sigmoid growth function was employed to model the dynamics of average basal area increments. Our study shows how the sigmoid growth model behaves with data derived from a broader growth period, i.e., the later stand development phase. The upper asymptote of the first and second derivatives of the model were chosen to measure the growth dynamics. We advocate using the basal area instead of the diameter (both at breast height) to estimate the growth dynamics. The simple explanation is that the same radial/diameter increment results in different net wood productions on trees with varying diameters in the initial period (assuming all other conditions are the same). This study was conducted to compare the growth dynamics of the vast majority of Polish larch provenances, i.e., the progeny of 21 populations growing at an altitude ranging from 50 to 650 m above sea level. The study's second goal was to compare the growth of the larch provenances with the Schober deterministic model (the stand volume and yield tables) [12] to find differences between the observed and expected growth dynamics. In many cases, the long-term experimental trials showed greater current stand growth dynamics compared to the yield table values [13]. The outcomes of this trial have severe implications for the legitimacy of using deterministic models for stand growth predictions.

Materials and Methods
The measurements were carried out in one of the five larch provenance experimental sites established in 1967 in cooperation with leading forestry science centres localized at the Siemianice forest experimental station (referred to as "Siemianice" hereafter). This study took measurements including and up to 2011. A detailed description of the seed collection, seedling production, and the localization of the maternal stands can be found in the study by Kocięcki [14]. The newest summary with detailed information about the maternal (seed) stands, the location and the basic climatic measures are provided in the study by Szymański and Wilczyński [7]. The experiment was primarily established following a completely randomized block design. Again, the details can be found in studies conducted by Rzeźnik [15,16], together with a comprehensive summary of the Polish research effort focused on larch populations [17,18]. Data from 86 plots of 21 provenances were chosen (see Table 1 for details), of which measurements of at least 17 inventories from 1966-2011 were available. The quadratic mean diameter was calculated: where dbh was the diameter at breast height, and n was the number of trees in the plot. The average basal area (g) was then calculated from the following formula: The calculation was repeated for each plot within each inventory. Next, the following four-parameter growth model was fitted with the function drm of the drc R package [19]): where b was the growth rate, c was the lower asymptote, d was the upper asymptote, e was the inflection point, x was the stand age (in years), and y was the average tree basal area at age x. Parameter c was fixed (c = 0) due to the known value of g at age 0. Thus, only three parameters had to be estimated. As a result, the four-parameter growth model was reduced to a three-parameter model. The mentioned function drm, employed for the parameter estimation, employs a general optimization based on the Nelder-Mead, quasi-Newton, and conjugate-gradient algorithms (function optim from the R package stats). Ritz et al. [19] referred to it in the description of the drm function and mentioned that for a continuous variable, such as in our study, it is reduced to a least square estimation. Figure 1a shows an example of a growth curve based on Formula 3. For each equation (one for each of the plots), the first and second derivatives were determined. To achieve this, we employed function D implemented in stats R package. The upper asymptote of the first derivate ( Figure 1b) indicates a point in time (X-axis) when the average basal area growth rate achieves its maximum. In turn, the upper asymptote of the second derivative ( Figure 1c) indicates a point in time (X-axis) where the basal area growth acceleration takes place [20]. The extrema mentioned above were determined by the function optimize (stats). We accomplished the same for the data from the deterministic model (the stand growth and yield tables [12]-based on the Schober tables [21]) for all three site indexes and the same range of development (stands up to 45 years old). The expected (average) growth curve was constructed from all the data points derived from the experiment. The coefficient of determination ρ 2 was calculated as follows: where, for each plot, j = 1, . . . , 86, the stand age i = 1, . . . , 46 g ij is the empirical basal area, g ij is the modelled value (Equation (1)) of plot j in age i, and g j is the empirical arithmetic mean for plot j.

Figure 1.
An example of the growth rate model (Equation (3)) (a) and its first (b) and second (c) derivatives. The sigmoid growth model (a) indicates cumulative growth for the investigated variable as a direct value. Here, the average basal area is a function of time. The first derivative (b) indicates the growth rate (physical speed) and achieves the maximum in the inflexion point of the growth model (Equation (3)). After that, the growth rate starts to decrease. The second derivative (c) is a speed of shift in the growth rate (physical acceleration). The extrema of that function indicate maximum acceleration (upper asymptote) and deceleration (lower asymptote).  (3)) (a) and its first (b) and second (c) derivatives. The sigmoid growth model (a) indicates cumulative growth for the investigated variable as a direct value. Here, the average basal area is a function of time. The first derivative (b) indicates the growth rate (physical speed) and achieves the maximum in the inflexion point of the growth model (Equation (3)). After that, the growth rate starts to decrease. The second derivative (c) is a speed of shift in the growth rate (physical acceleration). The extrema of that function indicate maximum acceleration (upper asymptote) and deceleration (lower asymptote). Since the experiment was established, windfalls and other random factors damaged many plots. Due to the violation of orthogonality (not all provenances were represented in each block), blocking factors could not be included. A one-way analysis of variance was used to compare the parameters and the extreme values of the derivatives from different provenances. All statistical analyses were performed using R 4.4.0 [22]. For the analyses of variance, we employed function aov from package stats. The analyzed variables were right-skewed, and the data were log-transformed to respect the assumption of normality. The residual distributions were compared with the normal distributions using the Shapiro-Wilk test using the function shapiro.test in the stats package. The homogeneity of the variance assumption was tested using the Levene test using the function leveneTest in the car package [23]. Post-hoc analysis was performed with the Tukey HSD test using the function HSD.test in the agricolae package [24]. The growth parameters of the experimental plots were compared with those estimated from the Schober deterministic model with a one-sample t-test. For all tests, an α = 0.05 significance level was assumed. The statistical analyses were performed in R 4.1.0 [22].

Results
At age 46, the average quadratic mean diameter ranged from 24.3 to 42.1 cm (31.1 on average ±3.49 standard deviation, SD). The sigmoid growth model used in our study (Equation (3)) explained 98% ± 2 SD of the observed variance. The supplementary materials describe the model performance among all analyzed plots in detail. The analysis did not reveal any significant differences between the basic growth model parameters, namely, the growth rate, upper asymptote, and the inflection point (Table 1) among the compared provenances ( Table 2). The upper asymptote, maximal growth rate, and maximal acceleration rate (Table 1, Figure 2) showed much higher variance (the coefficients of variation were 35%, 34%, and 40%, respectively) in comparison with the growth rate and inflection point (13% and 15%, respectively). A significant difference was found for the maximal acceleration rate. The Góra Chełmowa provenance achieved the highest mean value, but it was only significantly different from the under-performing provenances from Czerniejewo, Krościenko, and SzczytnaŚląska (Table 1, Figure 2). Table 3 reveals a difference in the sigmoid model estimates between all provenances taken together (the mean value for 86 plots) and those calculated for the deterministic model, i.e., the stand volume and yield tables (site index I). The latter had a significantly greater growth rate of 6% and a lower upper asymptote of 47%. When the extremes of the first and second derivatives are compared, i.e., the maximum growth rate and the maximum acceleration rate, the mean values for the provenances are higher by 77% and 69% from the deterministic model, respectively. The mentioned extremes occurred at similar stand ages (a difference of <1 year) in both the sigmoid and deterministic models, and there was no basis for rejecting the null hypothesis ( Table 3). The expected acceleration curve (Figure 2a), the growth rate curve (Figure 2b), and the vast majority of points representing the experimental plots, were higher than the curves based on the data from the deterministic model.

Growth Dynamics among Larch Provenances
The applied sigmoid growth model significantly explained the observed variability of the basal area growth dynamics. Only one growth dynamic parameter, i.e., the maximum growth acceleration, was indicated to significantly explain the variability among the tested larch provenance. The Post-hoc test showed that only the Góra Chełmowa provenance overpermed the three others (Czerniejewo, Krościenko, and SzczytnaŚląska) in this measurement. Szaban et al. [25] found that the Czerniejewo provenance had the highest average wood density. This means that the provenance indicated as an underperformer assimilated more carbon per volume unit to some extent. Moreover, wood samples from trees originating from the Czerniejewo provenance had the highest compressive and flexural strength [26,27] out of the six investigated provenances. Hence, when the aspects of wood's technical features are considered, growth dynamics might not be a good indicator of the overall provenance performance.
The results from an older experiment, established in 1949 in Rogów, are in line with our results. The stands representing different provenances reached the maximum rate of volume growth when they were 30 years old [17]. The Góra Chełmowa provenance outperformed the rest of the provenances for diameter from an early age [16]. The families of the population tested in Kórnik also showed more outstanding performance than yield table values [12], but only by 10% [28]. We also observed that the Grójec and Myślibórz provenances were above average. This is similar to Andrzejczyk and Bellon's study [29]. On the other hand, the Bliżyn and Skarżysko provenances belonged to a group of leaders as far as growth is concerned after the first 30 years of experiment in Rogów [29] but were still below average in Siemianice.
Based on the data from the IUFRO 1944 experiment, Giertych indicated that the fast growth in height was most often strongly related to low quality [30]. One of the over-performers (the Góra Chełmowa provenance) showed the most outstanding growth dynamics. Although previous reports from the studied experimental site revealed that it was characterized by poor stem forms from early developmental phases [16]. However, Rożkowski et al. [28] did not indicate significant differences in stem form among the nine families originating from the Góra Chełmowa population but mentioned that features were not compared among different provenances. In other words, the absence of significance among the families says nothing about this provenance if we focus on the interpopulation level. The Myślibórz provenance, which was also characterized by a high growth rate, was indicated to be the best in terms of the frequency of individuals with straight trunks [16,29]. Hence, the growth rate is not a good predictor for stem quality.

Expected Growth Dynamics
Generally, the volume and yield tables from Szymkiewicz [12,21] assumed a much lower growth dynamic than the data from the Siemianice experiment. The maximal growth rate was achieved at a similar time (inflection point, ca. 35-year-old stands) when the experimental trial and curves were based on the deterministic model (the stand volume and yield tables, site index I) were compared. Similarly, the maximum growth acceleration was achieved in similar ages when the data from our provenance trial and deterministic model were compared (ca. 22-year-old stands). Although, the average value of the maximal growth rate of the experimental stands, as well as the maximal growth acceleration, were almost twice as high and even three times higher for the Góra Chełmowa provenance. Szymkiewicz [12], who questioned the usefulness of Schober's deterministic model (the stand volume and yield tables) for Polish populations, was correct. By comparing the results published in the study by Kowalczyk and Matras [31] with the Schober deterministic model, a similar trend is noticed. They compared provenances to those planted in Siemianice, which grew under other conditions (a mountain mixed coniferous forest in the Sudetes). Ten Polish provenances growing in Lubawka achieved an average diameter at a breast height (DBH) of 45 years, ranging from 24.2 to 27.1 cm. The expected deterministic model value for site index I was 23.8 cm. This trend persisted until the last inventory documented in the study by Kowalczyk and Matras [31]. At 65 years, Polish provenances' average DBH ranged from 31.9 to 39.1 cm while the expected value was 31.2 cm.
Since we observed considerable differences between the deterministic model (based on historical data [21]) and the average growth rate of various larch provenances, two potential explanations may be derived. Firstly, the data used to develop the deterministic model came from populations characterized by low growth dynamics, less productive sites, or both simultaneously. Secondly, it is possible to state that climate change might impact the growth we observe currently. For example, an increased supply of carbon dioxide and an increase in nitrogen deposition could be one of the causes. The latter leads to nutrient disorders in trees [32]. A study on Pinus sylvestris L. indicated slight additional nitrogen deposition but this positively affected growth in the short term [33]. Simultaneously, in the cited study, the authors claimed that phosphorus availability was a limiting factor and that the unbalanced N:P ratio may determine growth in the future. The tripartite interaction of increasing atmospheric CO 2 concentrations, rising temperatures, and nutrition is a complex problem. From the study by Lukac et al. [34], we understand there is comprehensive knowledge of how singular factors affect trees. Still, we have little knowledge about the long-term influence of the interaction of multiple factors, and sometimes, there is no easy way to interpret the results from already established experiments. Given the projections about future species distributions modelled by Dyderski et al. [35], the larch will be the "loser" among species because the distribution of its optimal habitats will be restricted due to changes in abiotic conditions. Unfortunately, the above projection does not consider at least three crucial factors: species composition, genetic plasticity, and biotic factors. The greater biodiversity of tree species in stand composition increases stand productivity [36] and stability [37]. The latter is also dependent on population genetic plasticity. Dostálek et al. [38] indicated low to moderate genetic differentiation among the studied mountain populations across the native distribution area of L. decidua. In the Dostálek et al. study [38], larch populations from the Carpathian region were relatively genetically homogenous, but there were significant differences between the regions of the Alps. Here, we focused on the differences in the growth of different larch provenances (including those from the lowlands). We did not find any clear evidence for the growth rate dependency from provenances. Considering biotic factors, we understand that insects drive tree mortality [39]. Still, the impact of insect outbreaks in the future is supposed to be almost unpredictable [40].
Like other woody species, larch reacts to water deficits by decreasing wood production [41], which was also clearly shown in the study by Szeligowski [42]. Their dendroclimatological analyses of Polish larch provenances showed that all analyzed data (n = 14) indicated depression in tree ring sequences simultaneously. Szeligowski indicated, for some of the tested provenances, a slightly greater relative ring width during drought for Kowary, Grójec, Konstancjewo Góra, Konstancjewo Tomkowo, Rawa Mazowiecka, Bliżyn, Szczytná Sląska, and Skarżysko. From those provenances, our study revealed that the first three populations had a higher maximum growth rate than the average for all the plots. Hence, if the results shown in the Szeligowski study [43] can be generalized, three provenances could be recommended because of the resistance to drought and growth rate. Based on our results, this part of the discussion may be taken only as an attempt to explain what we observed. However, further study under the recommendations stated by Giertych [6] must be conducted, i.e., to compare the performances of larch populations among the gradient of site conditions.

Final Remarks
The results from the German provenance experiments showed that it is too early for a general inference about the growth of particular provenances [44]. The author mentioned above showed an example when the productivity of provenances underperforming in the first half of stand development became the most productive in the second half. A similar, much earlier conclusion was already provided by Giertych [6]. They indicated that a generalization based on experiments might be only preliminary even after 50 years of observations. We are convinced that the maximal growth acceleration has already been achieved (Figure 2a). By comparing with the deterministic model (stand volume and yield tables), we can state it is very likely that the maximum growth rate was also achieved in most of the experimental plots (Figure 2b). Nevertheless, it is currently difficult to predict the level of the upper asymptote of the average basal area. That is because it is too early to determine the point of the second derivative minimum and the maximum growth deceleration (Figure 1c). In other words, it is not reasonable to estimate this value when the observation range did not exclude this growth phase. This point is characteristic of the slow growth phase [45] when the stand matures. We provide the following example to emphasize a possible misinterpretation of the results. The growth curves for the deterministic model have upper asymptotes equal to 584, 411, and 281 (the I, II, and III site indices) when fitted to the data up to the age of 45. The exact estimates calculated with data up to the age of 140 are 2322, 1716, and 1219, respectively. It is definitively too early to estimate this value for the experimental stands.

Conclusions
A question arises concerning the provenance of the best growth characteristics. Among the 21 provenances evaluated, only one was statistically significantly superior (α = 0.05), i.e., the Góra Chełmowa provenance in the maximal growth acceleration. Many others were also above average. Although, for the other parameters (e.g., the growth rate, inflection point, maximum growth rate, and time of its occurrence), we have no justification for the claim about more significant among-provenance variance compared to the intra-provenance variance. In other words, we do not know which factors could explain the differences in the growth between populations, but it is probably not genetic. The present study only explored the growth dynamics. Future studies should focus on quality measures, which are still crucial in the wood industry. On the other hand, the focus of the research to find families or populations more resistant to drought or other unfavourable environmental factors, or simply plastic against the environmental variability, would be the solution from the point of view of projected climate changes.
The method used here could be applied to other long-term experiments. It is easy to implement and provides informative and straightforward stand growth characteristics. This paper has demonstrated a significant difference between the observed growth dynamics observed in the experimental trial of Polish provenances compared with the expected values from the deterministic model (the stand volume and yield tables). As a result, the growth assessments based on the deterministic model may underestimate the actual values for the conditions of Central Poland.