Determination of the Environmental Factors that A ﬀ ect the Growth and Survival of Greek Fir Seedlings

: Forests in the montane-Mediterranean zone have only recently began to be a ﬀ ected by wildﬁres, therefore the knowledge necessary for restoration projects is missing. The aim of the study was to assess the e ﬀ ects of factors related to seedling attributes, weather conditions and site suitability on seedling performance. The characterisation of sites was based on bedrock and soil clay content as well as pre-ﬁre vegetation. Apical growth and survival of seedlings was monitored for four years in Parnitha National Park. The parameters of a linear mixed model were estimated using annual apical growth of seedlings surviving in the end of the study as the dependent variable and type of site, rainfall, initial seedling height and age as explanatory ones. A quantile regression model using all the data available was estimated for each year of study, taking into account only initial height and site type as well as a logistic regression model of survival. The ﬁndings indicate that the growth of Greek ﬁr seedlings depends on May rainfall mediated by soil clay content, which in turn depends on bedrock, which is consistent with the “inverse texture hypothesis”. Sites with low soil clay content were always more beneﬁcial for survival, which was stronger a ﬀ ected by summer–autumn rainfall. In both contexts, drought stress due to soil clay content fades with increasing age. Sites that were not ﬁr dominated prior to ﬁre proved unsuitable also for planting ﬁr seedlings. A minor part of the observed variability could be associated with the initial height of seedlings, especially for seedlings showing high rates of apical growth.


Introduction
The natural vegetation in the Mediterranean area, as elsewhere, is characterised by a zonation defined by altitude. While low-and middle-altitude zones are mostly dominated by fire-adapted vegetation, this is not the case for the high-altitude zones. The Greek fir is a species endemic to south-central Greece that forms forests covering large areas in the montane-Mediterranean zone of its geographic distribution, typically between 900 and 1400 m asl, while isolated individuals are present from 600 to almost 2000 m asl [1]. This zone was hardly affected by wildfires in the past, but it has been increasingly affected in recent years. According to [2], 95,000 ha located higher than 900 m asl were affected by wildfires between 1983 and 2008; to the best of our knowledge, the first large scale wildfire affecting a Greek fir forest occurred in 2000 in Mainalo, Peloponnese.
The species is susceptible to fire and natural regeneration cannot lead to a recovery of the forest [3][4][5]. Fir is shade tolerant and regenerates poorly in open sites. Practically no natural regeneration was observed in the study site several years after the fire, even under the canopy, mortality was very high. The first two years are crucial for the survival of newly established plants; the chances significantly Land 2020, 9, 100 3 of 15 species. Nine 0.15 ha permanent plots were established in the burned area (Figure 1), each containing 55-75 seedlings planted in spring 2008; these were unevenly spaced due to the rugged terrain. The seedlings were 2 years old at the time of planting. For the purpose of the present study the plots were monitored from 2009 to 2012 and grouped in three types: flysch_fir: previously dominated by Greek fir, flysch substrate, clay content 23%-29%, altitude 970-1050 m asl (three plots), ls_fir: previously dominated by Greek fir, limestone substrate, clay content: 29%-43%, altitude 1020-1065 m asl (four plots), flysch_mixed: previously not dominated by Greek fir, flysch substrate, clay content 17%-18%, altitude 910-930 m asl (two plots). The values of the soil clay content were retrieved from [21]. Soils are mostly shallow and especially in limestone areas, rock outcrops are apparent. The sites varied in exposure, but this variable was not used in the grouping of sites and analysis since it did not improve model performance and explanatory power.
235 m asl, about 8 km SE of the study area ( Figure 1). The long term  average annual precipitation is approximately 430 mm and long term average temperature 16.5ºC. Precipitation amount in the study area is reported to be approximately twice as much as in Dekelia [22]. Monthly weather data for the study period were also obtained from this station. The amount of rainfall in the study area is closely related to the amount in the adjacent area, and thus, these values reflect the year-to-year variation in the study area as well. No interpolation was performed to account for altitude, since only data from this single weather station were suitable for the purpose of the analysis, and thus, the resulting values would only reflect the variation of the original ones. The most important rain events for seedling growth and survival are the ones that take place from late winter onward to late autumn, since they contribute to the building up of water reserves in the soil in spring until the cessation of growth in summer, and further provide the necessary water to overcome the typical Mediterranean summer dry period. Since our fieldwork showed that growth ceased in July (see also [23]), at the latest, we experimented with amounts of rainfall from February to July as an indication of water availability during the growth period [24]. Finally, the amount of May precipitation was used, since it resulted in the models with the best fit, based on the maximising R 2 and minimising AIC values. A similar approach was undertaken with respect to survival in order to identify the timing of rainfall producing the observed patterns. The rainfall during the critical months of the study period can be found in Table 1.  Long term climatic data for this area are available from the Dekelia meteorological station of the Hellenic National Meteorological Service, located at the southern foothills of the mountain at 235 m asl, about 8 km SE of the study area ( Figure 1). The long term  average annual precipitation is approximately 430 mm and long term average temperature 16.5ºC. Precipitation amount in the study area is reported to be approximately twice as much as in Dekelia [22]. Monthly weather data for the study period were also obtained from this station. The amount of rainfall in the study area is closely related to the amount in the adjacent area, and thus, these values reflect the year-to-year variation in the study area as well. No interpolation was performed to account for altitude, since only data from this single weather station were suitable for the purpose of the analysis, and thus, the resulting values would only reflect the variation of the original ones. The most important rain events for seedling growth and survival are the ones that take place from late winter onward to late autumn, since they contribute to the building up of water reserves in the soil in spring until the cessation of growth in summer, and further provide the necessary water to overcome the typical Mediterranean summer dry period. Since our fieldwork showed that growth ceased in July (see also [23]), at the latest, we experimented with amounts of rainfall from February to July as an indication of water availability during the growth period [24]. Finally, the amount of May precipitation was used, since it resulted in the models with the best fit, based on the maximising R 2 and minimising AIC values. A similar approach was undertaken with respect to survival in order to identify the timing of rainfall producing Land 2020, 9, 100 4 of 15 the observed patterns. The rainfall during the critical months of the study period can be found in Table 1. The plots were surveyed one or more times in summer and autumn of each year. Apical growth of each seedling was recorded at the mm level. Only the complete annual growth was taken into account in the statistical models of growth presented here. Two kinds of statistical models of growth were built: a. a linear mixed model (lmm) using apical growth as dependent variable and type of site, rainfall, initial seedling height and age as explanatory ones was fitted to the data of all years; the dataset used included only the individuals that were surviving in autumn of 2012 and b. a series of quantile regression models, one for each year of study, using all the data available for that year (i.e., including individuals that did not survive until the autumn of 2012). Since age and rainfall were identical within the same year, they were not included in the analyses. The same variables were used to build logistic regression models of seedling survival within each year. Similar to quantile regression growth models, only seedlings that had survived the previous year were taken into account and the response variable reflected the chance of survival until October of the year in question.
All the statistical analyses were performed in R environment [25]. We used lme4 [26,27] to perform the linear mixed effects analysis. Seedling height at planting and age (without interaction term) were treated as fixed effects. The effects of seedling age on growth were entered in the model to account for interannual differences in growth due to age, and thus, not to be confounded with differences in rainfall. Intercepts for individual seedlings as well as intercepts and slopes for the amount of rainfall by site type were modelled as random effects. The correlation between random slopes and intercepts was removed [28]. In sum, the following linear mixed model (lmm) was fitted to the data: , where: i: individual seedling, k: year, j: site type and growth ikj : apical growth in mm, height i : initial seedling height in mm, age ik : seedling age in years, Rainfall k : May rainfall recorded in Dekelia weather station in mm.
The profiled confidence intervals for the fixed effect parameter estimates and the variance parameters of the random effects were also estimated with the lme4 package. The values of apical growth were square root transformed to avoid heteroscedasticity that was detected in preliminary analyses run with the untransformed values. The BOBYQA optimiser was used throughout and the number of possible iterations the optimiser was allowed to have was set to 100,000. The residual plots were visually inspected without detecting deviations from homoscedasticity or normality. Additionally, the Brown-Forsythe Levene-type procedure of the lawstat package [29] was used to assess differences in residual variance among individuals and site types; in both cases no heteroscedasticity was detected. R 2 values [30] were estimated by means of the MuMIn package [31].
A quantile regression approach was adopted to analyse growth within single years [32]. The analyses were carried out with package quantreg [33]. The R 1 measure was estimated as a goodness of fit indicator. Note that, although its values lie between 0 and 1, R 1 is a local goodness of fit measure in this case, while e.g., Ordinary Least Squares R 2 is a global measure and these two are not comparable. Due to the estimation procedure R 1 values are always lower than OLS R 2 values that would be estimated from the same dataset. Finally, the logistic regression models were fitted by the stats package [25] and the goodness of fit was assessed by the Hosmer-Lemeshow test by means of the ResourceSelection package [34]. The number of groups was set to 10, since the number of samples was less than 1000 [35].

Results
The values of seedling height at planting followed a normal distribution (Kolmogorov-Smirnov p < 0.001) with mean 110.3 mm and standard deviation 37.6. The growth data of the seedlings surviving in 2012 (Table 2) were used to estimate a linear mixed model with height at planting, age, rainfall and site type as explanatory variables (Table 3).  Table 3. Estimated parameters of the lmm model of apical growth (mm) and associated confidence intervals and standard deviation and associated confidence intervals. Age: seedling age in years, initial height in mm, rainfall: May rainfall in mm. The effects of initial height and age on annual growth were positive and statistically different from 0; however, they contribute little in explaining the variation of seedling growth. R 2 m, marginal R 2 , which refers to the fixed effects, equals to 0.016. On the contrary, R 2 c, conditional R 2 , which refers to the whole model, amounts to 0.634, showing a strong effect of individual variation and rainfall, conditional on sites. The slopes reflect the change in square root apical growth (in mm), triggered by a 1 unit change in the quantitative explanatory variables. Most of the variance accounted for by the random effects is attributed to the intercepts reflecting individual variability. Secondly, it is attributed to intercepts of rainfall conditional on site types, although at the margin of being statistically different from 0, and less so to the different slopes of rainfall conditional on site type.

Fixed effects
The values of May rainfall were found preferable for estimating the model over a range of other month values or sums of month values (see Materials and Methods). The analysis reveals that the effects of rainfall on apical growth differ among sites. In the case of flysch_mixed, the slope and intercept are the lowest of all site types. With respect to flysch_fir and ls_fir, the intercept is larger in the former case and the slope is larger in the latter. This pattern suggests that increasing rainfall has a much stronger positive effect on seedlings planted in limestone sites than on seedlings planted in flysch sites, while during years of low rainfall flysch sites support higher growth. The estimated parameters point at the slopes crossing at slightly over 20 mm May rainfall recorded in the Dekelia weather station (which, as a rule of thumb, corresponds to twice as much in the study area): below this threshold growth is higher in the flysch sites, whereas the opposite is true if rainfall exceeds it ( Figure 2). This low rainfall condition occurred seven times in the years between 2000 and 2013, exactly in half of the occasions, according to the data provided by this weather station.
Land 2020, 9, x FOR PEER REVIEW 6 of 16 Table 3. Estimated parameters of the lmm model of apical growth (mm) and associated confidence intervals and standard deviation and associated confidence intervals. Age: seedling age in years, initial height in mm, rainfall: May rainfall in mm.  The apical growth per year of all the surveyed seedlings is summarised in Table 4. Seedlings that survived until the end of the study exhibited higher growth in all cases. The difference in apical growth between the seedlings that did survive until 2012, compared to those that did not, was The apical growth per year of all the surveyed seedlings is summarised in Table 4. Seedlings that survived until the end of the study exhibited higher growth in all cases. The difference in apical growth between the seedlings that did survive until 2012, compared to those that did not, was highly significant in 2009 and 2010 (Wilcoxon test, p < 0.001) and less so in 2011 (Wilcoxon test, p = 0.05). The parameter values of the estimated quantile regression models for each year's median growth (tau = 0.5), taking into account all seedlings, are presented in Table 5. Table 5. Parameters and confidence intervals for quantile regression models of apical growth (mm), tau = 0.5. Correlation coefficients for site groups concern differences with group ls_fir. The R 1 values estimated for the models of the first two years are very low (0.05 and 0.08 for the first and second year respectively), which rise in the last two years (0.  A more complete picture arises by looking at the graphs of the coefficients estimated for a more full range of quantiles ( Figure 3). While the slope values for initial height are generally close to zero, they depart therefrom in the last two years but only for the higher quantiles. Both the positive and negative effects of site type flysch_mixed are stronger in higher quantiles and this pattern is stronger in the last two years as well. The effect of site type flysch_fir exhibits the most complicated pattern. While the positive effect increases in higher quantiles in year 1 and 2, it shows signs of being non-monotonous thereafter, with the intermediate quantiles being negatively affected and the low and the high ones being closer to zero.  The parameters of the fitted logistic regression models of survival show a positive effect of initial seedling height, a usually non-significant difference between ls_fir and flysch_mixed and a positive effect of flysch_fir compared to ls_fir of varying magnitude ( Table 6). The strongest effects of site on survival were observed in 2011, a year of low summer rainfall. Table 6. Logistic regression parameter values (statistical significance ***: p<0.001 ***, **: p<0.01 , *: p<0.05) and Hosmer-Lemeshow goodness of fit test (g=10). Correlation coefficients for site groups concern differences with group ls_fir. estimate sd The parameters of the fitted logistic regression models of survival show a positive effect of initial seedling height, a usually non-significant difference between ls_fir and flysch_mixed and a positive effect of flysch_fir compared to ls_fir of varying magnitude ( Table 6). The strongest effects of site on survival were observed in 2011, a year of low summer rainfall. Table 6. Logistic regression parameter values (statistical significance ***: p < 0.001 ***, **: p < 0.01, *: p < 0.05) and Hosmer-Lemeshow goodness of fit test (g = 10). Correlation coefficients for site groups concern differences with group ls_fir. In order to systematically describe the effect of rainfall depending on the low or high clay content soils of the sites that were previously dominated by Greek fir, we regressed the estimated coefficient reflecting the difference between ls_fir and fysch_fir in the case of the median regressions of growth and the logistic regressions of survival (Tables 5 and 6) on the rainfall of each year and seedling age in the respective year. In this case, the response variable was not growth or survival per se, but the coefficients reflecting the difference in growth or survival between these two site types. The models best describing the relationships were:

Fixed
[Growth b] = −68.96 * log[age] − 0.38 * [May rain f all] + 111.07 (adjusted R 2 = 0.98), indicating a linearly declining effect of site on growth between ls_fir and flysch_fir with increasing May rainfall. At high rainfall values, the curve crosses the y = 0 threshold, marking that high soil clay limestone sites become more favourable for growth than low clay flysch ones. The logarithmic transformation of age shows that the declining difference (and the associated quicker threshold crossing) due to increasing age is gradually levelling off. However, the difference in the goodness of fit with the regression model estimated with the untransformed values is not great and while the general idea of declining difference with increasing age can be trusted, the exact form of the decline is only indicative.
[survival b] = −1.89 * [age] − 0.02 * [cumulative June-October rain f all] + 13.07 (adjusted R 2 = 0.95), indicating a linearly declining difference of the effects of site on survival odds between ls_fir and flysch_fir with increasing summer-early autumn rainfall and age. Due to the very low rainfall correlation coefficient compared to the intercept in this case, it is doubtful whether a threshold crossing, in the sense that ls_fir sites will become more favourable for survival than flysch_fir ones, will ever occur.

Discussion
The starting point of the present study was the hypothesis that the annual apical growth and the survival of outplanted Greek fir seedlings would be affected in the same way by factors that distinguish seedlings, sites and years among each other, namely seedling height at planting, geological substrate and soil attributes and prefire vegetation of the site as well as rainfall. The hypothesis was in most cases confirmed, although the observed patterns were less straightforward in the case of growth than in the case of survival.
The positive effect of individual attributes related to seedling vigour, represented in this study by the values of height at planting, was confirmed in both contexts (Tables 3, 5 and 6 and Figure 3) [36]. However, initial height contributes very little to explaining the variation in annual apical growth. Its importance rises in the last year of the study and the higher quantiles of apical growth, suggesting an increasing importance of individual vigour for growth once the seedlings become established. On the contrary the association of survival to initial seedling height becomes weaker at the later stages. The association of initial height to growth at later phases both in general [10,37] and under Mediterranean climatic conditions in particular [38] was often found to be positive. Seedling height is a common measure of seedling quality and evidence suggests that taller seedlings have an advantage that is retained during the later growth periods. However, this effect is often weak compared to the effects of environmental variables and sometimes even undetectable, especially in a rugged Mediterranean terrain [12].
The effects of site type on growth were partly congruent with the effects on survival. In the flysch_mixed site group, in which fir was subdominant prior to fire, growth rate was positively affected in the first year after planting compared to ls_fir and negatively thereafter (Table 5 and Figure 3). Probably resprouting evergreen and deciduous broadleaf plants had a nurse effect in the first year, which fades later on [39] (see also [36]) and/or the very low clay content of the soils in these sites drastically reduced water stress during the early phase of establishment compared to the other sites, from then on the unsuitability of the sites compared to the others surfaced [20]. The negative effect is stronger in higher quantiles indicating a genuine limiting factor from year 2 of the study onward. This interpretation is supported by the findings concerning survival ( Table 6). The odds of surviving in flysch_mixed were no different from ls_fir despite the unsuitability of the former group of sites. In the low rainfall summer of 2011, the odds of survival became significantly larger than in the more suitable, (in other respects) ls_fir, which should be attributed to the low soil clay content.
The effects concerning apical growth of site group flysch_fir compared to ls_fir were more variable, even showing an alternating pattern ( Table 3, Table 5 and Figure 2, Figure 3). Limestone sites were less beneficial for survival but in some years, they were found to be more beneficial for growth than sites on flysch. Different bedrocks lead to different soil nutrient environments, and also to different soil physical properties having a range of variable effects on outplanted seedlings [15]. In this area, soils developed from limestone parent material tend to be more fertile than those developed from flysch and contain a higher proportion of clay [40]. The effects of site type at later years are stronger in the seedlings showing moderate growth rates, at the low and high end the observed patterns are affected by different variables. Although effects on survival vary in strength among years, flysch_fir is always more beneficial than ls_fir ( Table 6).
The findings concerning the importance of rainfall timing on growth was consistent with other recent studies on fir populations (Greek fir and others) in Greece [17,41,42]. These studies pointed to May rainfall as the principal limiting factor for fir growth in southern Greece. The critical period shifted from May to early summer with increasing latitude. The fact that the strength of the effect of site on survival tends to decrease with increasing June to October rainfall is consistent with the well-known timing of drought stress in Mediterranean climate areas. The discrepancy between the timing of the rainfall that matters for growth and the timing that matters for survival suggests that high May rainfall might become a sort of ecological trap. Seedlings that have grown vigorously in late spring might be more vulnerable to an unusually dry summer.
Geological substrate has been found to exert a strong influence in postfire regeneration [43]. Postfire management should take into account fine scale differences in water related stress that may have serious implications on restoration [44]. While in this study the positive effects of increasing rainfall were confirmed in general, as anticipated, they differed among site types. The strength of the positive effect of the low clay flysch sites compared to high clay limestone ones is a declining function of age and summer rainfall. As the seedlings become established or summer rainfall increases, the odds tend to become more even among these site groups, partly in contrast to growth where the site effect is related to May rainfall. The combination of the findings concerning site and rainfall confirmed the inverse texture hypothesis in this case and shed light to the conflicting patterns observed. In 2011 and 2012, when May rainfall was above the estimated threshold, site group flysch_fir had a negative effect on growth compared to ls_fir, while the opposite is true in low rainfall 2010 (Tables 2 and 4). 2009 does not seem to fit in this picture, however it is the first year after planting and seedlings are expected to be harder affected by water availability due to limited root growth. The rainfall levels would have been sufficient to support higher growth in limestone sites in this year, should the seedlings be well established. In short, the unsuitability of site flysch_mixed is confirmed, while the initially positive effects of flysch_mixed on growth is reversed in later years; the evidence points towards mechanisms concerning water relationships and firm establishment with growing age. The shifting of importance hierarchies among management or site factors with increasing age is often confirmed in different settings [45]. On the other hand, while flysh_fir is constantly more favourable with respect to survival, although the size of the gap fluctuated depending on rainfall, it loses its edge with respect to growth in rainy years.
In water-limited environments, soil physical properties are often overlooked despite their importance for the partitioning of available water [46]. The idea that coarse-textured soils support higher plant growth compared to fine-textured ones under conditions of water stress is increasingly being supported by empirical evidence [47][48][49][50][51]. The main mechanism proposed to account for this effect is the deeper percolation of the limited rainwater leading to the preservation of higher amounts of water for transpiration, while in finer soils water is retained in the surface, thus resulting in a higher proportion of it being lost to evaporation [18]. A number of studies have shown that the diverging patterns of water partitioning between transpiration and evaporation generated by soil texture lead to the inverse texture effect [52]. Alternatively, this effect was attributed to the ability of clay to retain water, making it inaccessible to plants [47]. Additionally, unlike the hydraulic behaviour at or near saturation, hydraulic conductivity of fine textured soils is higher compared to coarse-textured ones at low soil water content levels, leading to a quicker drying out.
In areas characterised by high interannual rainfall variability, the advantage for growth can shift from the fine-textured soils in wet years to coarse-textured ones in dry [48]. Crossovers and thresholds concerning the effects of amount of rainfall on vegetation growth do not concern only or necessarily mean annual values, but may concern patterns of seasonality in precipitation specific for different plant functional types [49]. Although the inverse texture hypothesis was originally conceived for dry areas, it seems that it holds widely for plants facing drought stress; note that the objects of the study were seedlings and summer drought is a well-established bottleneck in Mediterranean areas, and also in high elevations [17]. The Mediterranean is known for its geographic variability including the occurrence of different soils within limited distance. The climate in the Mediterranean is also known for its interannual variability; the combination of both leads to the variable responses of plant growth among the years and sites observed in this study.
While the linear mixed model contributes to the understanding of the inter-and intra-annual variability in growth rates of seedlings surviving until the last year of the study, it still does not account for a large part of it ( Table 3). The quantile regression models (Table 5) for annual growth of all seedlings surviving in each year also account for a small part of the intra-annual variability. There are two sources of variability that can be considered responsible therefore. The high variability among seedlings leading to low within year explanatory power persists, even in later years when the number of surviving ones was a fraction of the number planted in the first place. The seed material used to produce the seedlings comes from a natural forest and Greek fir is known for its high intraspecies variability [53][54][55]. Thus, we hypothesised that variability in growth rates among seedlings would be high and the hypothesis was indeed confirmed. The lmm that allowed for a different intercept for each Land 2020, 9, 100 12 of 15 seedling fared better than quantile regression models and, indeed, a substantial standard deviation value was associated to the different intercepts of the individual seedlings. Allowing also for different slopes would have probably fared even better but the model became overparameterised and the algorithm could not converge. The importance of individual seedling traits, along with environmental variables, in models of outplanted seedling performance is observed also in settings dissimilar with ours [37]. The second potential source of variability not accounted for in our models is the intrasite variability. The rugged terrain includes a finer scale variability (at the scale of meters), including variability in soil conditions, than reflected in the broad site groupings adopted in this study, which was not accounted for [12,56].

Conclusions
Rainfall, either annual or within parts of the year or the growth period is commonly used as an index of water availability or severity of drought stress. The growth and survival patterns observed in outplanted Greek fir seedlings suggest that it is an appropriate index only for use as a between year comparison of the same or similar sites. If different sites are concerned, care should be taken to include soil variables in assessing water availability. To the best of our knowledge, this is the first time that the inverse texture hypothesis was considered in the context of reforestation as well as in the context of a drought sensitive phase of the plant life cycle in an environment otherwise not severely stressed by drought. The growth of Greek fir seedlings in Mt. Parnitha National Park was heavily dependent on May rainfall mediated by soil type. It was not rainfall or soil type but their interaction that lead to year-to-year variation. At the first stages of establishment and at years of low rainfall, firs grew better on flysch, and at later stages and rainy years they grew better on limestone. On the other hand, seedlings in flysch sites always had better chances of survival; also in this case, the effects faded with increasing age. More pronounced water stress results in higher mortality in limestone sites, but by being more fertile, these same sites provide opportunity for better growth for those that survive and become established unless in years of low rainfall. The crossover value of rainfall lies well within the usual range of rainfall in the study area, however at the low end of the usual precipitation values of the geographical distribution of fir dominated forests. Once the seedlings are well established, initial height dominates the growth pattern for quick-growing individuals and site characteristics dominate the pattern of moderately growing ones. Thus, the exact position of the rainfall threshold is a moving target since age has to be taken into account. A large part of the observed variability may be attributed to intraspecies variability not accounted for by height at planting or to small scale terrain variability.