The Effect of Environmental Factors on the Nutrition of European Beech (Fagus sylvatica L.) Varies with Defoliation

Despite being adapted to a wide range of environmental conditions, the vitality of European beech is expected to be significantly affected by the projected effects of climate change, which we attempted to assess with foliar nutrition and crown defoliation, as two different, yet interlinked vitality indicators. Based on 28 beech plots of the ICP Forests Level I network, we set out to investigate the nutritional status of beech in Croatia, the relation of its defoliation and nutrient status, and the effects of environmental factors on this relation. The results indicate a generally satisfactory nutrition of common beech in Croatia. Links between defoliation and nutrition of beech are not very direct or very prominent; differences were observed only in some years and on limited number of plots. However, the applied multinomial logistic regression models show that environmental factors affect the relationship between defoliation and nutrition, as climate and altitude influence the occurrence of differences in foliar nutrition between defoliation categories.


Introduction
Meteorological and climate conditions affect vegetation directly and indirectly. Direct effects include responses to temperature; indirect effects occur primarily as soil-mediated phenomena, such as the influence of precipitation on soil moisture regimes [1][2][3][4]. Probably the greatest threat to existing forest ecosystems in Europe comes from a changing climate [5]. Extreme climatic events such as droughts are thought to be important in initiating changes in forest ecosystems [6]. According to IPCC [7] predictions, extreme climate events such as heat waves and long-lasting summer droughts will occur more frequently. Such climate changes will negatively impact the stability, structure and biodiversity of forest ecosystems throughout Europe [8]. The region of southeast Europe is one of the most vulnerable regions in Europe when it comes to climate change impacts, primarily through intensified severity and duration of droughts and heat waves. As these impacts should be stronger and faster than on the rest of the continent [9,10] this region is an ideal model for studying the future impact of changing climatic conditions.
A prerequisite for healthy growth and balanced metabolism is the maintenance of adequate concentrations and relatively stable ratios of nutrients in plant tissues [11]. Concentrations of nutrients and their relationships in the leaves of forest trees are important indicators of their functioning. It is crucial to account for nutrient limitation when studying the forest response to climate change [12], as nutritional status has a significant and diverse impact on tree vitality [13]. The concentration of nutrients in leaves depends on many process, we investigated whether tree nutrition plays a role in the loss of vitality of beech trees (and vice versa); and if the link between tree nutrition and vitality loss varies from year to year due to interannual climate variations.

Foliar Nutrient Concentrations
Mean N foliar concentrations in the study area level were within normal range [40] in both defoliation categories and across sampling years. Repeated Measures ANOVA indicated that there was no difference in N concentrations between LD and HD trees, while there was a difference between individual years (F(2, 556) = 30.4, p < 0.0001). In 2018, significantly lower N concentrations were recorded for both defoliation categories compared to 2019 and 2020 (Table 1). Table 1. Results of Holm pairwise comparisons of foliar nutrient concentrations and sample dry mass within a defoliation category and between sampling years (NS-not significant, * p < 0.05, ** p < 0.01, *** p < 0.001).

Defoliation Category
Year Although there is no significant main effect of defoliation category on P concentrations, the pairwise comparisons between LD and HD trees indicated significant differences in 2019, when LD trees had a higher concentration than HD trees, but not in 2018 and 2020 ( Figure 1). Significant differences in P concentrations (F(2, 556) = 74.9, p < 0.001) were also observed between individual sampling years (Table 1) with a decrease during the study period. However, values of P on the plot level were mostly normal throughout the study period except for the year 2020, when P concentrations were slightly below normal range on 53% of the research plots ( Figure S3). No significant differences were established for potassium (K) between defoliation categories and mean concertation values stayed within the normal range in each year ( Figure 1).
Repeated Measures ANOVA showed no significant main effect of defoliation category on Mg concentrations. Nevertheless, pairwise comparisons between LD and HD trees revealed significant differences in 2018, with HD trees having higher Mg concentrations ( Figure 1). Mean Mg concentrations were within normal range in each year, although they were lower in 2019 and 2020 compared to 2018 for HD trees.
Similar to other elements, there was no significant main effect of defoliation category on calcium (Ca) concentrations (F(2, 556) = 3.4, p = 0.06). Pairwise comparisons, however, showed that HD trees in 2018 and 2019 had significantly higher Ca concentrations than LD trees ( Figure 1). For these years, values for both groups fall within normal Ca concentration range for beech. In 2020, both categories had significantly higher concentrations than in previous years, bordering on surplus (luxury) values, but no significant difference was found between them ( Figure 1).
Although there was no significant main effect of defoliation category on dry leaf mass, the pairwise comparisons between LD and HD trees indicated significant differences in 2018, when LD trees had a higher dry leaf mass than HD trees ( Figure 1). However, there was a statistically significant main effect of the sampling year on dry leaf mass (F(2, 556) = 56.3, p < 0.001) and pairwise comparisons showed differences between all years for both defoliation categories (Table 1).
On the plot level, mean Ca concentrations were mostly in the normal range during the first two years, while excessive Ca concentrations were observed on all plots located Mean foliar concentrations and sample dry mass (dots) with confidence intervals (vertical bars) for LD (green color) and HD trees (orange color) on all plots within sampling years. Dotted lines represent lower and upper critical value for normal range of foliar N concentrations for common beech according to Mellert and Goettlein [40]. Only significant differences in foliar concentrations between LD and HD categories determined by the Holm-Bonferroni method are indicated (* p < 0.05, ** p < 0.01, *** p < 0.001).
On the plot level, mean Ca concentrations were mostly in the normal range during the first two years, while excessive Ca concentrations were observed on all plots located in the alpine region in 2020 ( Figure S5), and most of the alpine region plots had insufficient P concentrations ( Figure S3). In the same year, excessive K concentrations were found on eleven plots (39%), most of them located in the continental region ( Figure S4). N values are rather plot dependent and relatively stable in time ( Figure S2), while Mg concentrations show a very broad range of values, independent of the area, defoliation category, or sampling year ( Figure S6) Overall, the share of plots with differences in foliar concentrations between defoliation categories was relatively low for most elements during the entire sampling period ( Table 2). The highest share of plots with differences in foliar concentrations between defoliation categories was determined for Ca (21.4%) in 2018, followed by Mg (17.9%) in the same year.

Foliar Nutrient Ratios
N/P ratios were especially high and above the upper limit for both LD and HD trees in 2020 at the study area level, mostly due to very low P values that approached deficiency values ( Figure 2). Repeated Measures ANOVA did not show a significant main effect of defoliation category on N/P ratio. However, pairwise comparison between LD and HD trees indicated significant differences in N/P ratios in 2019, with values out of optimal range for HD trees due to low P concentration values ( Figure 2). The main effect of the sampling year was significant (F(2556) = 13.4, p < 0.001) with an increase in N/P ratio during the monitoring period. Mean foliar ratios (dots) and confidence intervals (vertical bars) for LD (green color) and HD trees (orange color) on all plots within sampling years. Dotted lines represent lower and upper critical value for normal range of foliar ratios for common beech according to Mellert and Goettlein [40]. Only significant differences in foliar ratios between LD and HD categories determined by the Holm-Bonferroni method are indicated (* p < 0.05, ** p < 0.01, *** p < 0.001). Mean N/P, N/K and N/Ca ratios on the plot level were mostly in normal range during the entire study ( Figures S7-S9). However, in 2020, most of the plots located in the alpine biogeographic region had excessive N/P and insufficient N/Ca ratios.
The highest share of plots with differences between defoliation categories (5 out of 28, or 17.9%) was determined for N/Ca in 2018 and N/P in 2019. The share of plots in which we found differences hardly passed 10% for other ratios and years (Table 4). Mean foliar ratios (dots) and confidence intervals (vertical bars) for LD (green color) and HD trees (orange color) on all plots within sampling years. Dotted lines represent lower and upper critical value for normal range of foliar ratios for common beech according to Mellert and Goettlein [40]. Only significant differences in foliar ratios between LD and HD categories determined by the Holm-Bonferroni method are indicated (* p < 0.05, ** p < 0.01, *** p < 0.001). No significant differences were found for N/K ratio between defoliation groups or sampling years and values were within normal range in all years.
Repeated Measures ANOVA identified a significant main effect of defoliation categories on N/Ca ratios (F(2556) = 109.9, p < 0.0001). Pairwise comparisons confirmed differences in N/Ca ratios between defoliation groups in 2018 and 2019 ( Figure 2). These differences were dictated primarily by Ca concentrations, which showed higher values in HD trees in 2018 and 2019. Sampling year also had a significant effect (F(2556) = 13.4, p < 0.001) and pairwise comparison established that values within defoliation categories were lower in 2020 compared to 2018 and 2019 (Table 3).
Due to low Mg concentrations in 2019 and 2020, the N/Mg ratios were mostly in the surplus range in those years ( Figure 2). No significant effect of defoliation categories on N/Mg ratios was found. However, pairwise comparisons showed a significant difference between LD and HD trees in 2018, with higher values for LD trees. Table 3. Results of Holm pairwise comparisons of foliar nutrient ratios within a defoliation category and between sampling years (ns-not significant, * p < 0.05, *** p < 0.001). Mean N/P, N/K and N/Ca ratios on the plot level were mostly in normal range during the entire study ( Figures S7-S9). However, in 2020, most of the plots located in the alpine biogeographic region had excessive N/P and insufficient N/Ca ratios.
The highest share of plots with differences between defoliation categories (5 out of 28, or 17.9%) was determined for N/Ca in 2018 and N/P in 2019. The share of plots in which we found differences hardly passed 10% for other ratios and years (Table 4).

Effects of Environmental Factors on Differences in Foliar Nutrition
While the fitted MLR (multinomial logistic regression) models for N and P did not indicate a significant influence of any of the tested factors, we found that environmental factors influence the occurrence of differences in foliar concentrations of K, Ca, and Mg between defoliation categories. Altitude (alt) and mean annual temperature (MAT) are the only environmental factors that influence the occurrence of Differences in Foliar Concentrations (DFC) of K. Keeping all other variables constant, the odds for LD trees having higher K concentration is 0.5% higher for every meter increase in altitude. On the other hand, odds are 35% lower for HD trees to have higher K concentrations when the mean temperature increases.
Trees on higher altitudes are more likely to have higher Ca concentrations in HD trees (Table 5). An opposite effect can be seen with mean annual precipitation (MAP) where, keeping all other variables constant, with an increase in precipitation the odds are 90% higher for LD trees and 2% lower for HD trees to have higher Ca concentrations. With increasing maximum temperatures (Tmax) it is 1.3 times likely that HD trees will have significantly higher Ca concentrations.
Beech trees on higher altitudes are also more likely to have higher Mg concentrations in HD trees, but the effect on LD trees is not significant. Keeping all other variables constant, with a 1 • C increase in mean temperature it is 6.6 times likely that HD trees will have higher Mg concentrations, while at the same time the likelihood of LD having higher concentrations is negligible (Table 5). Table 5. Results of multinomial logistic regression for Differences in Foliar Concentrations (DFC) of potassium, calcium and magnesium and corresponding model AIC. Odds ratios are reported for each environmental factor while the confidence intervals are given in the brackets (* p < 0.1, ** p < 0.05, *** p < 0.01). The results of the fitted MLR models indicate that beech on higher altitudes is less likely to have significant differences in N/P ratios and increasing maximum temperature has the same effect (Table 6). A similar effect of maximum temperatures can be seen for differences in N/Ca. Holding all other variables constant, with an increase in precipitation the odds are 0.6% higher that HD trees will have higher N/P ratios and 0.8% higher N/Mg ratios. When mean temperatures and soil pH increase, it is 2.5 and 3.7 times, respectively, more likely that LD trees will have higher N/Mg ratios. Keeping all other variables constant, with a unit increase in PDSI, i.e., less drought, the odds are 83.5% higher that LD trees will have higher N/Ca ratios. The fitted MLR models for N/K did not indicate a significant influence of any of the tested environmental factors on the occurrence probability of differences in foliar nutrient ratios. Table 6. Results of multinomial logistic regression for Differences in Foliar Ratios (DFR) of potassium, calcium and magnesium. Odds ratios are reported for each environmental factor while the confidence intervals are given in the brackets (* p < 0.1, ** p < 0.05, *** p < 0.01).

Discussion
Assessing the nutritional status of trees using the values of element concentrations in needles or leaves is a common diagnostic practice in forestry [40][41][42]. In Croatia there have been only a few studies into the nutritional status of beech [41,43,44], however these studies were performed on a limited research area. The results obtained in this study indicate a generally satisfactory nutrition of common beech in Croatia. Mean concentrations on the study and plot levels were within the normal range for beech [40]. Similar results were reported by Ognjenović et al. [39] from a long-term monitoring based case study of beech foliar composition on one ICP Forests intensive monitoring plot.
Large deviations from the mean concentration values of longer periods in individual years have been reported for several mineral nutrients [45][46][47]. In contrast, foliar K concentrations were reported to be relatively stable [39], and this finding is repeated in our study, although concentrations of all investigated elements stayed mostly within limit values [40]. Very similar relationships of N and K to limit values has also been observed in other studies [48,49]. On the other hand, P foliar levels are a reason for concern at the European level [12,15]. While beech P nutrition on a study area level seemed to be sufficient [40], on average 33% of the plots had insufficient P nutrition during the study period ( Figure S3). Additionally, the ratio of nitrogen to phosphorus shows that beech nutrition in 2019 and 2020 was not balanced in this respect, which is partly caused by diminishing P, and partly by increasing N concentrations, respectively.
The fact that we found surplus N/Mg ratios in 2019 and 2020, both due to increasing N and decreasing Mg values, echoes results of other studies that report an unbalanced Mg nutrition on national [48,50] and European level [12].
Excessive Ca concentrations on an extremely high proportion of plots (90%) were determined in a study of beech nutrition in northern Spain, but since no symptoms of unbalanced nutrition have been identified in the field, the authors believe that the deviations of the obtained values may be based on ecological and geological differences between the studied forests and central European forests on the basis of which reference values were obtained [51]. This consideration is consistent with the hypothesis that the optimal leaf composition of each species depends on the specific ecological biogeochemical niche that the species occupies [52]. Taking into account the wide ecological amplitude of common beech, there is a possibility that the optimal range of concentrations of nutrients is adapted to the specific habitat in which a particular stand develops. Thus Sardans and Peñuelas [53] assume that plant species have a certain degree of flexibility in changing stoichiometry in response to changes in environmental conditions such as climate gradients.
Nutritional status and tree crown defoliation, as indicators of vitality, can provide a basis for monitoring the effects of long-term environmental changes on forest ecosystems [22,48]. According to Simon and Wild [54], if the concentration of a certain element remains in the normal range, a decrease in mineral nutrition should be regarded more as a consequence than the cause of damage. If, on the other hand, the concentrations are inadequate, we can suspect nutrition to be the cause of damage. Although we did not find a universal relationship between the defoliation categories and nutritional status of beech trees, subsequent comparisons revealed differences in the element concentrations and ratios for certain years. For instance, we found significantly higher P concentrations in LD trees in year 2019. On the contrary, a study in northern Spain recorded higher P concentrations in more damaged trees [51]. The most pronounced differences in concentrations between defoliation categories were determined for Ca, with HD trees having significantly higher concentrations in 2018 and 2019. Higher concentrations in HD trees were also recorded for Mg in 2018. Conversely, Jonard et al. [35] report that higher defoliation values are associated with lower Ca and Mg concentrations in a liming and P/K fertilization experiment of beech stands in Belgium.
A case study investigating interrelations of various common beech vitality indicators did not find significant links between foliar nutrient concentrations and defoliation [39]. It should be taken into account that (i) the study was performed on a limited area, (ii) element concentrations were mostly within the normal range, and (iii) the range of defoliation values was narrow and seldom exceeded 25%. The lack of association between N concentrations and defoliation values, also recorded in this study, Thimonier et al. [36] attribute to the same kind of limitations. Unlike Ouimet and Moore [55], who observed that an increase in K concentration in needles resulted in a decrease in defoliation in balsam fir stands, we found no significant differences in K concentrations between defoliation categories. Compared to fertilization experiments, the analysis of the nutritional status in natural stands is more demanding due to the high variability of data and the inability to control many factors that may affect the nutritional status. Ewald [56] states heterogeneous environmental conditions as the reason for the lack of significant relationship between nutrient concentrations and spruce defoliation status in a study conducted in Bavaria. On the other hand, with this study we do have a full picture of beech nutritional status in Croatia, including the insight into the significant variations between years and defoliation categories. Significant differences in concentrations only during certain years are not a phenomenon specific to this research alone. High year-to-year variability in nutrient concentrations is characteristic of nutritional studies [15,35,48,57].
The quantification of forest ecosystem response in a climate driven changing environment is fundamental for maintenance, enhancement and restoration of future forest ecosystem goods and services. The results presented in this study are a novel approach in investigating the effects of climate and environmental properties on the defoliation/nutrition relationship. Both tree nutrition and defoliation have been reported to depend on climate properties in various studies. A study by Braun [58] showed that air temperature and precipitation have a considerable impact on the foliar nutrient concentration. Jonard et al. [49] observed a positive relationship between precipitation and foliar Ca concentrations, while insufficient water supply during dry spells have a negative impact on the uptake of calcium, according to Bergmann [59]. Our results indicate that the increase in precipitation is better utilized in low defoliation trees since the odds are 90% higher for LD trees to have higher Ca concentrations with an increase in precipitation. Lukac et al. [60] state that elements with primarily biologically controlled cycles such as nitrogen can show a different reaction to temperature changes than elements whose cycles are controlled by both biological and geological processes (P and K) or predominantly by geological processes (Ca and Mg). Therefore, it is not surprising that we did not determine a significant effect of environmental factors on the occurrence of differences in N concentrations between defoliation categories. A study of beech nutrient dynamics along a precipitation gradient found that N concentrations remain constant with decreasing precipitation while P concentrations increase, resulting in a decrease in N/P ratio [61]. However, in this study increasing precipitation increases the odds of higher N/P ratios in high defoliation trees suggesting a disturbed uptake and utilization of P in beech with higher defoliation.
From our results it becomes quite clear that an unequivocal association of foliar nutrition and crown status should not be expected for beech. It seems plausible that various pressures influence this relationship, affecting the physiological status of beech trees in a multitude of ways. For instance, Gottardini et al. [22] found that various morphological and physiological characteristics of beech leaves and canopy such as leaf volume and photosynthetic activity, had a significant negative association with the extent of damage. Additionally, Ferretti et al. [62] showed that the 25%-defoliation threshold can be a reasonable approximation for tree classification indicated by the effect of slight and moderate variations in defoliation on tree growth, which is in contrast to the results obtained by Tallieu [63]. Highly differing research results on this topic are likely depending on the research area, management practices, climate influences and other unknown factors. Nevertheless, we established that environmental factors, especially climate, influence the occurrence of differences in foliar nutrition between defoliation categories.

Study Area and Plot Design
The ICP Forests Level I monitoring plots in Croatia are established on intersections of a 16 × 16 km grid that contain forest cover. These plots do not have a fixed area; rather, 24 trees are chosen for defoliation assessments using a cross-cluster system with six trees in each cluster [64]. Only plots with a minimum of five beech trees in the sample in 2018 were selected to ensure that beech was significantly represented in the mixture of tree species, resulting in total of 28 plots (Figure 3).

Study Area and Plot Design
The ICP Forests Level I monitoring plots in Croatia are established on intersections of a 16 × 16 km grid that contain forest cover. These plots do not have a fixed area; rather, 24 trees are chosen for defoliation assessments using a cross-cluster system with six trees in each cluster [64]. Only plots with a minimum of five beech trees in the sample in 2018 were selected to ensure that beech was significantly represented in the mixture of tree species, resulting in total of 28 plots (Figure 3).

Climate Data
Climate monitoring stations are generally situated at considerable distances from the forest research plots. Therefore, the data they provide are not always representative of the research locations. To overcome this, we used gridded data produced by regression kriging (RK), which is a hybrid method of interpolation carried out in four steps [68,69]. The method was validated with leave-one-out cross-validation, while the root mean square error (RMSE) was calculated between observed and interpolated values. Mean RMSEs are for mean monthly temperature from 0.5 °C to 0.9 °C, for minimum temperature from 1.1 °C to 1.5 °C, for maximum temperature from 0.7 °C to 1.1 °C, and for precipitation from 18 to 30 mm, averaged by months.
Monthly temperature and precipitation values from the gridded dataset on 1 km spatial resolution for Croatia [30] were used to calculate Mean Annual Temperature (MAT), mean annual minimum (Tmin) and mean annual maximum (Tmax) temperature, Mean Annual Precipitation (MAP) as well as the Palmer Drought Severity Index (PDSI) [70] and Standardized Precipitation Evapotranspiration Index (SPEI) [71]. Lower values of scPDSI

Climate Data
Climate monitoring stations are generally situated at considerable distances from the forest research plots. Therefore, the data they provide are not always representative of the research locations. To overcome this, we used gridded data produced by regression kriging (RK), which is a hybrid method of interpolation carried out in four steps [68,69]. The method was validated with leave-one-out cross-validation, while the root mean square error (RMSE) was calculated between observed and interpolated values. Mean RMSEs are for mean monthly temperature from 0.5 • C to 0.9 • C, for minimum temperature from 1.1 • C to 1.5 • C, for maximum temperature from 0.7 • C to 1.1 • C, and for precipitation from 18 to 30 mm, averaged by months.
Monthly temperature and precipitation values from the gridded dataset on 1 km spatial resolution for Croatia [30] were used to calculate Mean Annual Temperature (MAT), mean annual minimum (Tmin) and mean annual maximum (Tmax) temperature, Mean Annual Precipitation (MAP) as well as the Palmer Drought Severity Index (PDSI) [70] and Standardized Precipitation Evapotranspiration Index (SPEI) [71]. Lower values of scPDSI and SPEI indicate a stronger drought intensity while higher values indicate a higher degree of humidity. SPEI was calculated on a time scale of 3 months.

Tree Selection Procedures, Foliar Sampling and Analysis
Sampling was conducted annually from 2018 to 2020 during late July and August. At a distance of less than 50 m from the centre of the plot, a stratified sampling of trees based on the percentage of defoliation was carried out, forming two groups of five trees: (i) "LOW DEFOLIATION" (LD) trees with defoliation lower than 25% and (ii) 'HIGH DEFOLIATION' (HD) trees with defoliation higher than 25%. Defoliation assessments were performed during July and August according to the ICP Forests methodology [24] in steps of 5%, ranging from 0 to 100%. An absolute reference tree, defined as the best/healthiest possible beech tree, was used in the assessments, regardless of local factors that may affect the crown condition. Only dominant or codominant trees were selected, without any wounds or fruiting bodies of fungi on the trunk. In cases when trees changed their defoliation percentage and thus changed their group status (LD/HD), new, replacement trees were selected instead, and the discarded trees were not used in further sampling. Stratified sampling allowed us to circumvent the irregular distribution of tree defoliation status on plots and focus on determining whether concentrations/contents of nutrients in leaves differ in trees of different defoliation categories on each plot. Leaves were sampled from the upper third of the crown using a shotgun rifle or a rope climbing technique [72,73]. Samples from each tree were combined into a composite sample according to fresh weight and analyzed in the Laboratory for Physical and Chemical Testing (LFKI) of the Croatian Forest Research Institute. Upon arrival to the laboratory, samples of plant material were dried at 105 • C to constant mass, ground in a Fritsch Pulverisette 14 mill and prepared for analysis in a Milestone Ethos One microwave oven [74]. Concentrations of total nitrogen were determined on a Leco CNS 2000 elemental analyzer (LECO, 2000). Phosphorus (P) concentration was determined colorimetrically on a LaboMed UV/VIS spectrophotometer [74], and potassium, calcium and magnesium by atomic absorption on the Perkin-Elmer Aanalyst 700 absorption spectrophotometer [74].

Soil Sampling and Analysis
Soil sampling was performed in 2019 on five points located within each of the research plots. One point was located within each of the four clusters of trees that are assessed for defoliation during regular monitoring activities, and an additional fifth point was located in the centre of each plot. Soil samples were taken with a pedological drill from a depth of 0-10 cm, 10-20 cm, 20-40 cm, and 40-80 cm. Collected samples were pooled by sampling depth. Soil mechanical properties were determined by [75] and chemical parameters were analysed according to the following protocols and methods: • soil pH (CaCl 2 ) [76] • exchangeable K, Ca and Mg according to [77] • exchangeable acidity (free H + ) according to [78] • total N with an elemental analyzer Leco CNS 2000 [79] • available P and K with the AL method [80]; P by spectrophotometry using molybdate blue method, on UV/VIS spectrophotometer LaboMed and K directly from filtrate on flame photometer Buck scientific PFP−7 [81].

Data Analysis
Descriptive statistics of foliar concentration and ratios were performed for each element. Foliar data distribution was inspected visually and by Shapiro-Wilk test. Slight deviations from normality were determined, especially for potassium, calcium and magnesium. However, studies have shown that slight deviations from normality do not affect the rate of false positive results [82].
Therefore, a two-factor repeated measures analysis of variance (ANOVA) was used to determine the difference in foliar concentrations and ratios between defoliation categories (LD/HD) of common beech on a study area level. The data of all elements met the assumption of homoscedasticity, which was confirmed by Levene's test. The Holm-Bonferroni method [83] was used to perform pairwise comparisons to determine differences in foliar concentrations (i) between LD and HD categories in single years and (iii) across defoliation categories between sampling years.
To determine differences in foliar concentrations and ratios between defoliation categories on a plot level, we used the independent samples t-test [84]. The results were then categorized for each plot and sampling year to reflect three cases: (i) HD_higher if HD trees had significantly higher foliar concentrations (p < 0.05), (ii) LD_higher if LD trees Multinomial logistic regression (MLR) was used to determine which environmental factors influence the occurrence probability of a certain DFC/DFR category, except no_dif which was set as the reference category. Separate odds ratios, which are the exponentiation of model coefficients, were determined for all independent variables for each DFC/DFR category. The odds ratio of a coefficient indicates the occurrence probability of either HD or LD trees having significantly higher foliar concentrations over the probability of no difference between defoliation categories (NO_diff ). An odds ratio > 1 indicates that the occurrence probability of a certain DFC/DFR category relative to the occurrence of NO_diff increases as the independent variable increases, whereas an odds ratio < 1 indicates that the occurrence of a certain DFC/DFR category decreases. Of the potential independent variables, those with the highest variable importance were selected with the random forest algorithm [85]. The final model selection process was based on diagnostic diagrams and a procedure defined by Johnson et al. [86]. All analyses were conducted in an R programming environment [87].

Conclusions
Our final view is that links between defoliation and nutrition of beech in Croatia are not very direct or very prominent; this can be seen from (I) the fact that differences in nutrition of LD and HD trees were found for only a few nutrients (II) these effects are mostly not universal, but present only in some years and on a limited number of plots and (III) nutrient concentrations mostly stay within normal values regardless of tree defoliation status. However, it is clear that environmental factors, especially climate properties, do affect the relationship between defoliation and nutrition. Discovering mechanisms by which environmental factors affect foliar nutrition should complement current research on the defoliation-nutrition relationships.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants12010168/s1, Figure S1. Comparison of study area climate properties during the vegetation period with the climate normal (1980-2011) of the vegetation period (March -September); Figure S2. Foliar nitrogen concentrations on research plots; Figure S3. Foliar phosphorus concentrations on research plots; Figure S4. Foliar potassium concentrations on research plots; Figure S5. Foliar calcium concentrations on research plots; Figure S6. Foliar magnesium concentrations on research plots; Figure S7. Foliar N/P ratios on research plots; Figure S8. Foliar N/Ca ratios on research plots; Figure S9. Foliar N/K ratios on research plots; Figure S10. Foliar N/Mg ratios on research plots; Figure S11. Defoliation of common beech on research plots; Figure S12. Defoliation of common beech on all plots during the study period.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to legal reasons.