Variability and Relationship Between Phenological and Morphological Traits in Early and Late Pedunculate Oak

Abstract

Variability and relationship between phenological traits and leaf, acorn, and tree size morphometric traits were examined in early and late bud-flushing groups of the pedunculate oak population in the vicinity of Sremska Mitrovica, Serbia. According to the obtained three-year results, there were no significant differences in tree size and leaf morphometric characteristics between the early and late group. The effect of trees within phenological groups was statistically significant and considerable, especially in leaf blade width (lbw) and leaf area (la). Acorn length (acl) and acorn index (acinx) were significantly influenced by phenological group, achieving moderate contribution of phenological group to the total variation. There was a clear effect of phenological group on variation of examined phenological traits. The effect of year of monitoring on the same traits was not significant, but there was a clear effect of interaction between phenological group and year, especially in case of traits describing the period from bud swelling until the emergence of wrinkled leaves (f12a), unfolded leaves (f12b), and fully developed leaves (f12g) and their ratio with the period from 8 March until the emergence of wrinkled leaves: f12b/f02a and f12g/f02a ratio. Earlier phenology in warmer years is clearer in the early bud-flushing group than in the late one. Periods between different bud swelling and f12a, f12b, and f12g phenological phases were significantly shorter in 2017 (with a warm April) than in 2015 (with moderate temperatures in March and April) in the early group, but significantly longer in the late group. Examined traits were classified in six groups based on their factorial loadings with the first six principal components rotated by Varimax method, revealing strict distinction between traits by their original nature. In that sense, all examined groups of traits could be considered as informative in variability studies of pedunculate oak. The tree size traits (tree height and diameter at breast height) formed the separate, fifth group, suggesting no close relationship of these traits with any other examined characteristic. Both cluster analysis and PCA suggest distinct classification by trees’ phenology, but also considerable differences by the second principal component which is closely related to leaf size characteristics. The research should be continued on variability between populations and progenies, especially with respect to phenological and acorn morphometric traits. Understanding the phenological variations between early and late oaks could be essential for designing robust forest adaptation management strategies.

1. Introduction

Pedunculate oak (Quercus robur L.) is a noble broadleaved tree species and one of the dominant edificators in riparian forest communities of Southeast Europe. In Serbia, pedunculate oak is an autochthonous species whose distribution is conditioned by high ground water level and occasional flooding. The greatest stands of pedunculate oak are in Vojvodina province, in the riparian zones of rivers Sava and Danube and their tributaries, as well as in the central part of Serbia in the riparian zone of river Morava [1,2]. Those pedunculate oak forests are the most endangered forest complex by climate change in Serbia, since they are under the negative impact of rising temperatures as well as declining groundwater level [3,4].

The most preserved stands of pedunculate oak (29,000 ha) in Serbia are in the area named “Ravni Srem”, where several hundred years’ old trees were found in the nature reserves of “Smogva” and “Vratična” [5].

The broad areal and considerable variability of soil and hydrological conditions of habitats where pedunculate oak grow has resulted in distinct variability of morphological and physiological traits and numerous forms and varieties within the species [6,7].

However, the management of the pedunculate oak stands faces numerous difficulties in the last few decades due to the dieback of solitary and group trees at almost the whole area where this species is present [1]. Climate changes, drought in summer, insufficient precipitation and its unfavorable distribution, and unfavorable changes in water regime are the main causes of these problems. Considering the high ecological and economic value of pedunculate oak stands in Serbia, it is necessary to study more deeply and understand the ecology of pedunculate oak stands in order to develop effective restoration and conservation strategies for their survival.

In that sense, phenological research [8,9,10], as well as morphological [10,11,12,13] and ecophysiological studies [14,15,16], provide useful information on complex genetic diversity in nowadays considerably changed ecological conditions.

Genetic variability is a key factor in the adaptation of species to changes in environmental conditions and ongoing climate change. Thus, the low variability of certain species endangers their survival [17].

Among the most important adaptive traits, directly affected by climate change is phenology [18]. In the second half of the 20th century, the trends of temperature increase in forested areas of temperate climate were related to the earlier beginning and prolongation of the growing season [19]. According to Erdeši [20], considerable intraspecies variability of pedunculate oak includes phenological characteristics of early and late bud flushing. Based on phenological variability, two basic varieties of pedunculate oak in Serbia are defined: “early” (Quercus robur var. praecox Čern.) and “late” (Quercus robur var. tardiflora Čern.). The period between bud-flush phase between these two varieties is 2–5 weeks. These varieties also differ in growth characteristics, wood quality, and morphological traits [20]. Knowing individual phenological characteristics of host pedunculate oak plants can be effectively used in mitigating the damage caused by early-spring defoliators [21,22].

Understanding the dependence of phenology on abiotic factors is also important. Phenological events are primarily based on air temperature, although the importance of photoperiod, winter chilling, and early and late frosts is not yet known [19].

The set of a phenophase is under considerable genetic control, and it is useful in classification of basic material in individual selection regarding tolerance, growth, and fructification [23,24]. In that sense, Andrić et al. [25] suggest the establishment of clonal seed orchards from the late phenophorm pedunculate oak individuals, to insure the production of quality reproductive material.

The interaction between growing degree day and plant phenology is often present in phenological studies based on the influence of air temperature accumulation on plant developmental processes [26,27,28]. However, the studies on the effect of growing degree day on phenological traits and morphological traits are not common, especially for pedunculate oak.

Selection of pedunculate oak phenoforms better pre-adapted to the warming climate can be a feasible measure from ecological, genetic, and economic points of view; however, this can reduce genetic diversity and limit the maintenance of a population’s resilience and stability [29].

Studying leaf morphology can not only elucidate the phenological variability of pedunculate oak, but also provide useful information about genetic differentiation at intra- and interpopulation levels and can be a good basis for determination of species and lower taxons, as well as hybrids. Similarities between individuals within and between populations can provide insight into their historical relationship and common ancestry [30].

Leaves are plant organs where the processes of photosynthesis and transpiration take place. Variability of basic characteristics of leaves, such as shape, size, margin, petiole traits, leaf area, etc., can reveal the importance of the common influence of genetic factors and local environmental conditions [31]. Fu et al. [32] showed that leaf traits can vary depending on the climate, suggesting that they can express different phenotypes as the response to environmental variability, i.e., phenotypic plasticity [33]. Several studies conducted at the level of species showed that this plasticity is controlled by both genetic and environmental components [34,35].

Acorn traits are also used in population studies of oaks. Many studies in oak species have been performed regarding the variability of acorn size traits linked to geographical and aridity range [36,37]. Also, the relationship between acorn size traits, as well as between them and seedling growth and physiological characteristics, has been investigated by many authors [37,38,39,40,41,42]. Thus, these studies stress the importance of acorn size characteristics not only for population studies but also from the aspect of forest management.

We hypothesize that differences between and within phenophorm classes significantly affect their phenological pattern and leaf and acorn morphological traits, as well as differences between years, and that there is a significant relationship between these traits. In that sense, the aim of this work is to define the phenological patterns of early and late bud-burst groups of pedunculate oak test trees in three successive years, as well as to study variability of leaf and acorn traits in order to analyze variability of examined pedunculate oak trees from the area of “Gornji Srem”, Serbia.

2. Materials and Methods

For the study of variability and relationship of phenological and morphometric traits in pedunculate oak, thirty-six trees of population “Vinična-Žeravinac-Puk” near Sava River in the vicinity of Sremska Mitrovica, Serbia (44°55′50.36″ N 19°08′41.22″ E, altitude 89 m.a.s.l.) were selected. Two groups of trees were selected regarding their bud-burst phenology: seventeen trees with early bud burst (Q. robur var. praecox Čern.) (trees number 11, 13, 21, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, and 38), and nineteen of them with late bud burst phenology (Quercus robur var. tardiflora Čern) (trees number 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 15, 16, 17, 18, 23, 36, and 37).

Two tree size traits were measured: breast height diameter (DBH) [cm] and tree height (h) [m]; then leaf morphometric characteristics: leaf blade length (lbl) [cm], leaf petiole length (lpl) [cm], total leaf length (ll) [cm], leaf blade width (lbw) [cm], leaf area (la) [cm²], lbw/lbl ratio (lbw/lbl), and lpl/lbw ratio (lpl/lbw); then acorn morphometric traits: acorn length (acl) [cm], acorn width (acw) [cm], acorn index (acinx), and acorn volume (acve) [cm³]. The tree size traits were measured in 2014, while leaf and acorn morphometric traits were measured in the autumn of 2014, 2015, and 2016.

The basic features of the relief on which this population extends are flat terrain with alternating beams and depressions that are mostly parallel to the course of the Sava River. The geological base in the area consists of alluvial deposits of sand of different structure [43]. To describe the soil conditions in this study, pedological profiles were excavated. These populations are distributed on hydromorphic, semigley soils. Soil type is fluviatile meadow soil (humofluvisol) of heavier mechanical composition. Humofluvisols in the study area have a structure of the A-AC-C profile and they are characterized by a powerful humus horizon with a depth of up to 50 cm. Thickness of fresh leaf mass on the soil surface is about 2 cm. The transformation of organic matter is favorable. The percentage of humus is highest in the humus horizon and decreases with depth, ranging from 0.6% to 10.2%, and according to the classification, these soils are very humus-rich. On this area, the flood season lasts less than 30 days, and groundwater level lies at a depth of 160–180 cm.

The “Vinična-Žeravinac-Puk” population is on ecologically and economically optimal habitats for the growth and development of pedunculate oak stands. These are mixed stands of pedunculate oak with European hornbeam (Carpinus betulus L.) and narrow-leaved ash (Fraxinus angustifolia Vahl.) (Carpino-Fraxino—Quercetum roboris caricetosum remota) at the altitude of 82–87 m. Selected trees are in adult phase 40 to 80 years old, 50 m away from each other. The breast height diameter was calculated from two cross measurements of diameter at breast height performed by electronic caliper DP II (Haglöf Sweden 88221 Långsele, Sweden) with precision of 0.1 mm. The tree height was measured by VERTEX III V 1.5 and Transponder T3 (Haglöf Sweden 88221 Långsele, Sweden) digital instrument with the precision of 1 cm.

For the measurement of leaf morphometric traits, thirty leaves per tree were collected in August 2016 from short branches from the outer part of the lower half of the tree crown. Thirty healthy, undamaged leaves per tree were collected, herbarized, and prepared for further measurements.

Four leaf morphometric traits (lbl, lpl, ll, lbw) were measured according to reference [44], by ruler, with the precision of 1 mm, while three characteristics (la, lbw/lbl and lpl/lbw) were additionally included.

For measurement of leaf area, leaf samples were scanned, while the leaf area (la) was determined by computer program ImageJ ver. 1.49p [45].

For the measurement of acorn morphometric traits, fifty acorns per tree were collected (between mid-September and October 2014). The length (acl) and width (acw) were measured by digital caliper with the precision of 0.1 cm. The acorn index was calculated by formula: acinx = acl/acw; meanwhile, acorn volume was calculated by formula for ellipsoid as it was proposed by Aizen and Patterson [36]: acve = 4/3 π acl acw².

Monitoring of the bud-flushing phenology was performed on all 36 selected trees from 8 March (last date when all plants were dormant) until the end of May in three consecutive years (2015, 2016, and 2017), two times a week by the same person. Average monthly temperatures in January, February, March, and April were as follows: 2.4, 2.7, 7.1, and 11.9 °C in 2015; 0.9, 7.2, 7.6, and 13.8 °C in 2016; and −5.5, 3.8, 9.7, and 11.6 °C in 2017, respectively [46].

The determination of phenological phase was performed according to reference [30] where

• 1—bud swelling
• 1a—bud elongation
• 2a—appearance of the first wrinkled leaves 0.5–3 cm long
• 2b—leaves 3–5 cm long
• 2v—leaves 5–7 cm long, intense leafing
• 2g—leaves 7–10 cm long, optimum leafing.

Based on phenology monitoring data, the following traits were derived: f01—period from 1 January till phase 1 [day]; f02b—period from 1 January till phase 2b [day]; f12b—period from phase 1 till phase 2b [day]; f02g—period 1 January till phase 2g [day]; f12g—period from phase 1 till phase 2g [day]; f02a—period from 1 January till phase 2a [day]; f12a—period from phase 1 till phase 2a [day]; f12b/f02a—f12b/f02a ratio; f12g/f02a—f12g/f02a ratio. The starting date in bud-flushing phenology monitoring was 1 January [26].

Further, the growing degree day (GDD) was calculated by the following formula:

$$GDD={\sum }_n^l(T_m-T_{base})$$

where GDD stands for growing degree day, n for the first day of measuring period from 1 January, l for the last day of measuring period from 1 January, T_m for the mean daily temperature, T_base for the base daily temperature, where values for days that had mean daily temperature lower than 0 °C are not included in the formula (Pellis et al.) [26]. In this way, the following characteristics were derived: f01t: GDD (growing degree day) for the period from 1 January till phase 1 [°C]; f02at: GDD for period from 1 January till phase 2a [°C]; f12at: GDD for period from phase 1 till phase 2a [°C]; f02bt: GDD for period 1 January till phase 2b [°C]; f12bt: GDD for period from phase 1 till phase 2b [°C]; f02gt: GDD for period from 1 January till phase 2g [°C]; f12gt: GDD for period from phase 1 till phase 2g [°C]; f12at/f12bt—f12at/f12bt ratio; and f12at/f12gt—f12at/f12gt ratio.

Statistical Analysis

The examined traits were regarded as quantitative traits with normal distribution of frequences. Because tree size characteristics were measured only in one year, their variability is described by Box and Whisker plots. For leaf, acorn, and phenological traits, the variability was presented by interval of variation, F test from analysis of variance, coefficients of variation, and contribution to expected variance. The examined trees were separated into two phenological groups: early and late bud-flushing groups. For leaf and acorn traits, the two-way hierarchical analysis was performed, with Phenological group and Trees within Phenological group as controlled sources of variance For phenological traits, the two-way factorial analysis of variance was performed, with Phenological group, Year and Phenological group × Year as controlled sources of variation. The significance of differences between treatments was tested by Tukey’s HSD (Honestly Significant Difference) test at the significance level α = 0.05.

Principal component analysis was calculated for total means for all examined traits, and it was based on correlation matrix. The first three principal components were used for description of the relationship between trees and the factorial loadings of original variables (examined traits) with selected principal components rotated with Varimax method for grouping of examined traits.

Hierarchical cluster analysis with unweighted pair group method with arithmetic mean (UPGMA). The method was used for the analysis of agglomeration of examined trees according to dissimilarity between them based on Euclidian distance calculated from average values of trees, standardized by standard deviation of tree average values. The distance at which the clusters were formed was determined by scree test, based on agglomeration schedule. Statistical procedures were carried out with the STATISTICA 14.0.0 software package (TIBCO Software Inc., Santa Clara, CA, USA) [47].

3. Results

Based on one-way analysis of variance, there were no significant differences between examined phenological groups in breast height diameter (p = 0.233) and tree height (p = 0.166) (Figure 1).

According to the results of analysis of variance, there were no significant effects of differences between phenological groups on leaf and acorn morphometric traits, except for acorn length (acl) and index (acinx). However, the effect of trees within phenological groups was significant for all morphological traits. The coefficient of variation for phenological groups was generally weak for all morphological traits, where the highest value was recorded for acorn index (acinx). Higher were coefficients of variation for trees within phenological groups, especially for the length of leaf petiole (lpl) and length of leaf petiole/leaf blade width ratio (lpl/lbw), that were both 34.6%. According to the Tukey’s HSD test, the early bud-flushing group had longer leaves than the late group. The early group acorn was longer and narrower than the one for the late group, so its acorn index was higher than the index of the late group (

Table 1. Variability of examined leaf and acorn morphological traits of selected trees of Quercus robur based on two-way hierarchical ANOVA.

Traits (a)	Interval of Variation	F-Test				Coefficient of Variation (%)			Tukey’s HSD Test (c)
	Min–Max	FP (b)	pP	FT(P)	pT(P)	CVP	CVT(P)	CVres	Early Group		Late Group
lbl	7.880–14.713	1.188	0.284	4.815	0.000	1.634	13.852	35.421	10.890	a	10.274	b
lpl	0.387–1.968	0.146	0.705	6.800	0.000	0.000	34.640	71.838	0.683	a	0.717	a
ll	8.267–15.467	0.956	0.335	5.064	0.000	0.000	13.762	34.094	11.573	a	10.991	b
lbw	5.317–9.225	0.090	0.766	10.803	0.000	0.000	15.219	24.276	6.836	a	6.724	a
la	26.195–70.815	0.115	0.737	10.906	0.000	0.000	27.495	43.633	43.353	a	41.920	a
lbw/lbl	0.500–0.754	1.713	0.200	2.271	0.000	1.945	7.169	31.760	0.640	a	0.668	a
lpl/lbw	0.069–0.331	0.133	0.718	5.232	0.000	0.000	34.601	84.028	0.105	a	0.110	a
acl	2.144–3.464	7.772	0.009	99.619	0.000	5.698	9.228	6.571	3.023	a	2.773	b
acw	1.086–1.808	1.183	0.284	86.194	0.000	1.061	10.446	8.002	1.350	b	1.403	a
acinx	1.646–2.747	13.038	0.001	102.481	0.000	8.626	10.480	7.356	2.263	a	1.994	b
acve	1.513–5.489	0.003	0.959	88.312	0.000	0.000	26.529	20.076	2.968	a	2.954	a

^(a) Labels of traits: lbl—leaf blade length [cm]; lpl—leaf petiole length [cm]; ll—total leaf length [cm]; lbw—leaf blade width [cm]; la—leaf area [cm²]; lbw/lbl—lbw/lbl ratio; lpl/lbw—lpl/lbw ratio; acl—acorn length [cm]; acw—acorn width [cm]; acinx—acorn index; acve—acorn volume [cm³]. ^(b) Subscript labels in F (F-value), p (probability of F-test), and CV (coefficient of variation) value labels regarding source of variation: P—Phenological group, T(P)—Trees within phenological groups, res—Residual. ^(c) Difference between values with the same letter is not significant for α = 0.05.

).

The contribution of phenological group to the total expected variance clearly confirmed that there was no effect of this source of variation to the variation of most morphometric traits, except for acl and acinx; for those, the contribution was 20.0% and 31.2%, respectively. The contribution of trees within phenological groups was relatively weak in leaf morphometric traits (4.8%–28.4%), but rather high in acorn traits (46.1%–63.6%). The difference is due to high variation of leaves within trees (71.5%–94.8%), while the contribution of residual in acorn traits was relatively moderate (22.7%–36.7%) (Figure 2).

The two-way factorial analysis of variance was performed for Phenological group and Year as the main effects for examined phenological traits.

In this case, only traits that describe the period and growing degree days between 1 January and date of phenological phases 2a, 2b, or 2g were under significant main effect of Phenological group. The main effect of Year was not significant for any of the examined phenological traits, but the interaction Phenological group × Year was significant for all of them. The highest coefficient of variation for factor Phenological group was found for GDD traits, f12at and f12bt (50.0% and 40.3%, respectively), whilst the highest coefficient of variation for interaction Phenological group × Year was for traits that described the period from phase f1 to phase f2a (f12a) or phase f2b (f12b) (43.7% and 35.0%).

According to the Tukey’s HSD test, all phenological traits had higher values in late than in early bud-flushing trees in all three years, with few exceptions in that the difference was not significant. Also, the results of Tukey’s test additionally confirmed differences in reaction of examined phenological groups on examined years, as it was suggested by significant effect of Phenological group × Year interaction. For most phenological traits, the difference between 2016 and 2017 was not significant in the early group, or their values were significantly lower in 2017 than in 2016, as they were in f12a and several GDD traits such as f01t, f02bt, f02gt, f02at, and f12a/f12bt. However, in late group, phenological traits in general were significantly higher in 2017 than in 2016, except for f01, f02bt, f02gt, f02at, and f12a/f12bt ratio. Also, values of phenological traits were usually higher or not significantly different than in 2016 and 2017, but in late group, traits f12g, f12g, f12a, f12a/f12g, f12bt, f12gt, f12at, and f12a/f12gt were significantly higher in 2017, whilst in early group, f12a/f12b, f01t, f02bt, and f02at were significantly higher in 2016 than in 2015 (

Table 2. Variability of examined phenological traits of selected trees of Quercus robur for early and late phenological groups and examined years based on two-way factorial ANOVA.

Traits (a)	Interval of Variation		F-Test (b)						Coefficient of Variation (%) (b)
	Min–Max		FP	pP	FY	PY	FP×Y	pP×Y	CVP	CVY	CVP×Y	CVres
f01	78.941–98.263		11.430	0.077	2.405	0.294	17.098	0.000	0.082	0.037	0.043	0.045
f02b	93.588–122.526		20.689	0.045	1.281	0.438	25.913	0.000	0.130	0.019	0.050	0.042
f12b	13.588–36.737		2.364	0.264	0.247	0.802	60.076	0.000	0.238	0.000	0.350	0.193
f02g	110.235–136.737		35.665	0.027	2.863	0.259	10.319	0.000	0.105	0.030	0.029	0.041
f12g	31.294–50.947		1.802	0.312	0.317	0.759	36.674	0.000	0.100	0.000	0.191	0.136
f02a	87.569–115.474		21.198	0.044	1.261	0.442	21.947	0.000	0.131	0.018	0.049	0.046
f12a	7.569–27.842		2.388	0.262	0.189	0.841	60.357	0.000	0.300	0.000	0.437	0.241
f12a/f12b	0.531–0.755		1.209	0.386	0.123	0.891	33.935	0.000	0.038	0.000	0.140	0.104
f12a/f12g	0.226–0.545		1.728	0.319	0.003	0.997	51.960	0.000	0.157	0.000	0.316	0.188
f01t	513.330–776.111		16.298	0.056	8.992	0.100	10.202	0.000	0.137	0.121	0.057	0.080
f02bt	742.822–1193.153		50.004	0.019	1.707	0.370	8.868	0.000	0.238	0.035	0.056	0.084
f12bt	229.492–558.115		14.487	0.063	0.446	0.692	14.737	0.000	0.403	0.000	0.183	0.210
f02gt	1001.835–1466.824		42.082	0.023	1.475	0.404	7.943	0.001	0.181	0.024	0.046	0.074
f12gt	488.505–833.023		8.226	0.103	0.321	0.757	17.436	0.000	0.211	0.000	0.132	0.138
f02at	635.845–1060.465		47.533	0.020	2.261	0.307	8.988	0.000	0.243	0.049	0.058	0.087
f12at	122.515–410.372		11.372	0.078	0.333	0.750	17.657	0.000	0.500	0.000	0.261	0.271
f12a/f12bt	0.509–0.733		4.012	0.183	0.141	0.876	11.829	0.000	0.107	0.000	0.102	0.131
f12a/f12gt	0.243–0.493		13.280	0.068	0.177	0.850	8.068	0.001	0.306	0.000	0.141	0.225
Traits	Tukey’s HSD Test (c)
	Early2015 (d)		Early2016		Early2017		Late2015		Late2016		Late2017
f01	82.000	cd	78.941	d	80.000	d	98.263	a	88.158	b	85.789	bc
f02b	102.824	c	96.882	d	93.588	d	121.842	a	110.421	b	122.526	a
f12b	20.824	bc	17.941	cd	13.588	d	23.579	b	22.263	bc	36.737	a
f02g	122.177	c	110.235	d	112.529	d	135.526	a	128.842	b	136.737	a
f12g	40.176	b	31.294	d	32.529	cd	37.263	bc	40.684	b	50.947	a
f02a	95.353	c	91.706	cd	87.569	d	115.474	a	103.342	b	113.632	a
f12a	13.353	c	12.765	c	7.569	d	17.211	b	15.184	bc	27.842	a
f12a/f12b	0.643	c	0.718	ab	0.531	d	0.729	ab	0.669	bc	0.755	a
f12a/f12g	0.336	c	0.409	bc	0.226	d	0.457	b	0.373	c	0.545	a
f01t	544.653	d	676.293	b	513.330	d	737.240	a	776.111	a	607.453	c
f02bt	807.826	c	938.044	b	742.822	c	1162.926	a	1193.153	a	1165.568	a
f12bt	263.173	c	261.751	c	229.492	c	425.686	b	417.042	b	558.115	a
f02gt	1159.684	b	1191.745	b	1001.835	c	1453.020	a	1466.824	a	1440.476	a
f12gt	615.031	c	515.452	d	488.505	d	715.780	b	690.713	bc	833.023	a
f02at	710.999	c	833.096	b	635.845	c	1030.259	a	1060.465	a	1017.825	a
f12at	166.347	c	156.803	c	122.515	c	293.019	b	284.354	b	410.372	a
f12a/f12bt	0.639	b	0.607	b	0.509	c	0.686	ab	0.673	ab	0.733	a
f12a/f12gt	0.275	c	0.305	c	0.243	c	0.403	b	0.413	b	0.493	a

^(a) Labels of traits: f01—period from 1 January till phase 1 [day]; f02a—period from 1 January till phase 2a [day]; f12a—period from phase 1 till phase 2a [day]; f02b—period 1 January till phase 2b [day]; f12b—period from phase 1 till phase 2b [day]; f02g—period from 1 January till phase 2g [day]; f12g—period from phase 1 till phase 2g [day]; f12a/f12b—f12a/f12b ratio; f12a/f12g—f12a/f12g ratio; f01t: GDD (growing degree day) for the period from 1 January till phase 1 [°C]; f02at: GDD for period from 1 January till phase 2a [°C]; f12at: GDD for period from phase 1 till phase 2a [°C]; f02bt: GDD for period 1 January till phase 2b [°C]; f12bt: GDD for period from phase 1 till phase 2b [°C]; f02gt: GDD for period from 1 January till phase 2g [°C]; f12gt: GDD for period from phase 1 till phase 2g [°C]; f12at/f12bt—f12at/f12bt ratio; f12at/f12gt—f12at/f12gt ratio. ^(b) Subscript labels in F (F-value), p (probability of F-test), and CV (coefficient of variation) value labels regarding source of variation: P—Phenological group, Y—Year, P × Y—interaction Phenological group × Year, res—Residual. ^(c) Difference between values with the same letter is not significant for α = 0.05. ^(d) Labels for treatments of interaction Phenological group × Year: “Early” stands for the early bud-flush group, “Late” for the late bud-flush group; and the number stands for the year of monitoring.

).

According to the contribution of examined sources of variance based on two-way factorial ANOVA, Phenological group factor achieved dominant effect in f01, f02a, f02b, and f02g, while in f12a, f12b, f12g, f12b/f02a, and f12g/f02a, the dominant effect achieved Phenological group × Year interaction. The dominant contribution of controlled sources of variation in GDD traits achieved factor Phenological group. The contribution of the main effect of Year was rather poor in examined phenological traits except in f01t, while contribution of Residual effect ranged from 8.1 in f02b to 44.1% in f12a/f12bt (Figure 3 and

Table 3. Broad-sense heritability for leaf and acorn morphological traits and phenological traits.

Leaf Morphological Traits				Acorn Morphological Traits				Phenological Traits
Early Group		Late Group		Early Group		Late Group		Early Group		Late Group
lbl	0.329	lbl	0.071	acl	0.568	acl	0.726	f01	0.134	f01	0.000
lpl	0.438	lpl	0.150	acw	0.609	acw	0.647	f02b	0.013	f02b	0.000
ll	0.337	ll	0.075	acinx	0.643	acinx	0.698	f12b	0.000	f12b	0.000
lbw	0.278	lbw	0.284	acve	0.576	acve	0.677	f02g	0.000	f02g	0.226
la	0.311	la	0.262					f12g	0.000	f12g	0.000
lbw/lbl	0.041	lbw/lbl	0.151					f02a	0.061	f02a	0.014
lpl/lbw	0.085	lpl/lbw	0.167					f12a	0.000	f12a	0.000
								f12a/f12b	0.000	f12a/f12b	0.141
								f12a/f12g	0.000	f12a/f12g	0.000
								f01t	0.000	f01t	0.000
								f02bt	0.000	f02bt	0.444
								f12bt	0.156	f12bt	0.000
								f02gt	0.000	f02gt	0.508
								f12gt	0.000	f12gt	0.010
								f02at	0.000	f02at	0.486
								f12at	0.000	f12at	0.000
								f12a/f12bt	0.000	f12a/f12bt	0.183
								f12a/f12gt	0.000	f12a/f12gt	0.022

).

3.1. Principal Component Analysis (PCA)

The principal component analysis was based on the matrix of correlation coefficients calculated from the means of trees. According to the results of principal component analysis, the first three unrotated principal components described 70.2% of total variance. The scores of the first and third unrotated principal components (closely related to phenological and acorn traits, respectively) were used to form PCA biplot graph to relate examined trees (Figure 4). On this graph, the relatively clear distinction between early bud-flushing trees and the late ones along the first principal component is presented. Also, it seems that both early and late bud flushing trees are distributed along the third unrotated principal component, suggesting considerable contribution of acorn traits to the total variability.

To analyze the relationship between examined traits, the first six principal components were rotated by Varimax method. This rotation of principal components results in maximized variance of loadings within principal components. The selection of the first five principal components is done according to Keiser’s role, by which the principal components with eigenvalue higher than 1 are selected for further analysis. However, the sixth principal component was also included because the trait lbw/lbl had its highest loading with it (−0.739), while the eigenvalue of its principal component was close to 1 (0.966). Because this trait had low loadings with the other principal components and this principal component had low loadings with other examined traits, it seemed that the trait lbw/lbl carries information about differences between trees that is poorly described by other traits and could be expected to be grouped separately from other traits. The results of Varimax rotation were used for the analysis of the relationship between examined traits. Because there is no correlation between principal components, all traits that have their highest loading with the same principal component were in close relationship and grouped in the same group. All phenological traits are included in the first group of traits. They have their highest loadings with the first principal component. This principal component contributes the most to the total variance (53.3%), which suggests that examined trees differ dominantly by phenological traits. The only member of this group that is not a phenological trait is acorn index (acinx), but its loading with the first principal component is relatively low (−0.513). The traits of the second group are exclusively leaf morphometric traits: lbl, ll, lbw, and la, traits that describe the size of leaf, while the third group consists of acorn traits: acl, acv, and acve. The fourth group consists of traits based on petiole length: lpl and lpl/lbw, and the fifth group consists of tree size traits: d13 and h. The highest loading with the sixth principal component has only lbw/lbl ratio, describing the shape of the blade. It seems that no trait was found to be closely related to tree size traits, except for acinx that had relatively moderate loading with the fifth principal component (

Table 4. Factor loadings between original variables and the first six principal components rotated with Varimax method. The highest factorial loadings of original variables are underlined.

Original Variable (a)	PC1	PC2	PC3	PC4	PC5	PC6
d13	−0.243	−0.007	−0.125	0.156	−0.834	−0.008
h	0.218	0.031	0.029	0.144	−0.773	0.011
acl	−0.458	−0.343	0.570	−0.230	−0.348	0.124
acw	0.135	0.067	0.957	0.066	0.176	0.066
acinx	−0.513	−0.346	−0.359	−0.266	−0.481	0.053
acve	−0.059	−0.063	0.984	−0.045	0.013	0.101
lbl	−0.140	−0.949	−0.033	−0.016	−0.021	−0.211
lpl	0.128	−0.232	−0.004	0.902	−0.152	0.002
ll	−0.114	−0.949	−0.033	0.129	−0.044	−0.203
lbw	−0.068	−0.913	0.076	−0.041	0.017	0.361
la	−0.030	−0.931	0.065	0.020	0.020	0.185
lbw/lbl	0.158	−0.050	0.181	−0.089	−0.023	0.933
lpl/lbw	0.169	0.130	−0.027	0.916	−0.140	−0.112
f01	0.925	0.074	0.050	−0.007	0.028	−0.018
f02b	0.990	0.072	0.028	0.018	0.012	0.010
f12b	0.957	0.062	0.002	0.042	−0.005	0.038
f02g	0.979	0.078	−0.040	−0.033	0.023	0.000
f12g	0.890	0.070	−0.136	−0.058	0.014	0.021
f02a	0.987	0.073	0.023	0.066	0.025	0.011
f12a	0.963	0.065	−0.011	0.143	0.021	0.044
f12a/f12b	0.775	−0.024	−0.009	0.312	0.072	0.180
f12a/f12g	0.911	0.040	0.059	0.240	0.032	0.107
f01t	0.915	0.096	0.022	−0.013	0.038	−0.026
f02bt	0.989	0.073	0.023	0.029	0.012	0.009
f12bt	0.978	0.054	0.023	0.053	−0.005	0.030
f02gt	0.975	0.092	−0.038	−0.023	0.018	−0.020
f12gt	0.947	0.084	−0.072	−0.027	0.004	−0.015
f02at	0.986	0.079	0.015	0.087	0.035	0.005
f12at	0.976	0.062	0.009	0.154	0.030	0.027
f12a/f12bt	0.764	0.004	0.014	0.309	0.100	0.111
f12a/f12gt	0.919	0.038	0.081	0.240	0.039	0.104
Explained variance	16.514	3.904	2.443	2.230	1.748	1.241
Contribution to the total variance	0.533	0.126	0.079	0.072	0.056	0.040

^(a) Labels of traits: d13—diameter at the breast height (1.30 m); h—tree height; f01—period from 1 January till phase 1; f02a—period from 1 January till phase 2a; f12a—period from phase 1 till phase 2a; f02b—period from 1 January till phase 2b; f12b—period from phase 1 till phase 2b; f02g—period from 1 January till phase 2g; f12g—period from phase 1 till phase 2g; f02a/f12b—f02a/f12b ratio; f2a/f12g—f2a/f12g ratio; f01t—GDD (growing degree day) for the period from 1 January till phase 1; f02at—GDD for period from 1 January till phase 2a; f12at—GDD for period from phase 1 till phase 2a; f02bt—GDD for period 1 January till phase 2b; f12bt—GDD for period from phase 1 till phase 2b; f02gt—GDD for period from 1 January till phase 2g; f12gt—GDD for period from phase 1 till phase 2g; f12at/f12bt—f12at/f12bt ratio; f12at/f12gt—f12at/f12gt ratio; lbl—leaf blade length; lpl—leaf petiole length; ll—total leaf length; lbw—leaf blade width; la—leaf area; lbw/lbl—lbw/lbl ratio; lpl/lbw—lpl/lbw ratio; acl—acorn length; acw—acorn width; acinx—acorn index; acve—acorn volume.

).

3.2. Cluster Analysis

Cluster analysis was performed to relate examined trees. According to the scree test, five clusters were defined according to agglomeration step at the 6.1 of Euclidian distance. All early bud flushing trees were agglomerated in the second cluster. The only late bud flushing member of Cluster 2 was tree number 36, because it had similar leaf morphometric features as early trees, particularly high values for lpl, lbw, la, and lpl/lbw. All other trees are late bud flushing trees, agglomerated in Clusters 1, 3, 4, and five. Most of the late trees were classified in Cluster 5, tree number 5 is in Cluster 1, tree number 8 in Cluster 3, and tree 2 in Cluster 4. Thus, it seems that the main traits that defined clusters were phenological traits (Figure 5).

4. Discussion

4.1. Variability of Examined Traits

Variability of population and the entire species represents the basis of their stability [48]. In that sense, variability studies are crucial in protection and restoration of populations and species, in this case in studied Quercus robur population in Sava basin.

Four groups of traits were examined in this study: from tree size traits, then morphometric traits of leaf and acorn, to phenological traits. While morphometric traits, especially those of leaf, are widely used in pedunculate oak quantitative genetic and population studies, research on variability of phenological traits is rare. Because of that, the variability studies presented in this work could be considered as the basis for further studies on variability and the relationship of pedunculate oak traits.

Although our results showed no significant differences between the early and late phenological groups, it seems that early trees tend to grow greater diameter, but late trees tend to be higher. Greater timber is highly valued, but narrow rings that could be found in trees with smaller diameter could be interesting in carpentry. Nevertheless, the effect of phenological groups on tree growth should be the subject of further studies.

The examined leaf morphometric traits are shown to be poorly influenced by differences between examined phenological groups, although it seems that leaves in early trees were longer. The effect of Trees within the phenological group was statistically significant. Traits based on petiole length (lpl and lpl/lbw), dominated by high Tree within Phenological group coefficients of variation, but also had the highest Residual coefficients of variation, contrary to the other leaf traits. In that sense, the more reliable characteristic was the contribution to the expected total variance, by that the highest contribution of leaf blade width (lbw) and leaf area (la)—close to 30%. Also, broad sense heritability for lbw and la was found to be relatively moderate both in early and late group, while there was clearly higher heritability in early than in late group in lbl, lpl, and ll. Considering these results, lbw and la could be proposed as leaf traits that provide the most information on variability within the examined pedunculate oak population.

Examining the intra- and interpopulation variability of some pedunculate oak morphological leaf properties in 44 natural populations of Bosnia and Herzegovina, Ballian et al. [13,49] have determined statistically significant variability in each population for most of the studied leaf properties. These results are in accordance with the research of Franjić [50] on 20 populations of pedunculate oak from different bioclimatic areas of Posavina and Podravina in Croatia. The analysis revealed significant differences between populations, and differences between trees in the population, with the fact that the differences between trees in the population are greater than the differences between populations. The length of the leaf petiole was found to be the most variable leaf morphometric trait. Similar to findings presented Batos et al. [10] in five pedunculate oaks populations in Belgrade and Northern Serbia region, the characteristics with the highest variability are as follows: length of petiole, ratio between length of petiole to total lamina length, leaf area with petiole, and dry leaf mass. The leaf area had a strong positive correlation with maximum width of lamina, lamina perimeter with petiole, and length of lamina. In this research, results suggest that all of the studied leaf characteristics showed highly significant statistical differences between the trees, and also most of the characteristics showed differences that are significant for population differentiation.

Among acorn traits, variability of acorn length (acl) and acorn index (acinx) were significantly influenced by Phenological group, while the effect of Trees within phenological group was significant in all acorn traits. Acorns in early trees were longer but thinner, and the acorn index was higher in the early trees. Coefficient of variation of Phenological group was poor in all acorn traits, and the highest coefficient of variation of Trees within Phenological group and Residual were the highest in acorn volume (acve). However, the contribution of Phenological group was moderate for acl and acinx, but poor for acw and acve, while the effect of Trees within phenological groups was considerable for all acorn traits. This fact, with moderate to high heritability both in early and late trees, proposes these traits for further variability studies in pedunculate oak.

Woziwoda et al. [39] determined differences in dimensions and weight of North American red oak (Q. rubra) acorns. A weak correlation was found between the length and width of acorns in relation to volume and mass, which are strongly correlated. Popović et al. [38] determined that in species with large seeds, as in the case of Northern red oak, the seedling growth in the first growing season is closely related to seed size or the amount of reserve nutrients that are stored in the seed. The same correlation was obtained when it comes to root collar diameter. Devetaković et al. [40] found better performance of large than small pedunculate oak acorns in Northern Serbia, as they had higher percent of germination, greater dry weight of roots, and lower shoot to root ratios. Roth et al. [41], found that pedunculate oak (Q. robur L.) and sessile oak (Q. petrea L.) acorn size has a positive effect on germination percent, shoot growth, seedlings height, and survival rate in juvenile stages of development. Similar results are confirmed by Matić et al. [51], as well as Roth et al. [41,42].Considering their findings in positive correlation between acorn size and geographical range of the oaks in the eastern part of North America, Aizen and Patterson [36] suggested that acorn size can considerably influence the range size changes due to climate change. According to Bonito et al. [37], larger differences between Q. ilex populations in the growth characteristics they observed during the first growing season could be justified by the significant differences among acorn size. However, this effect of acorn size on growing and physiological seedling traits was not found in following years.

In temperate deciduous forests, the growth period can be defined as the period between bud burst in spring and leaf fall in autumn. Winter cooling and short photoperiod are known to delay early bud burst and prevent frost damage, while rising temperatures promote bud burst Reyer et al. [52]. Jovanović et al. [53] noted that favorable weather conditions at the beginning of growing period tend to push forward the beginning of some phenophases, but do not shorten but prolong them, especially in the case of the bud-flush phenophase. Generally, tree species populations from colder regions show the tendency to leaf earlier than those from warm regions. In that sense, southern oak populations tend to leaf earlier and suffer more damage from spring frosts, while those from regions with shorter growing season that finish their growth earlier avoid damage from autumn frosts [30,54].

Among several deciduous tree species, Vitasse et al. [55] described the phenological characteristics of oaks to be most flexible to changes in temperature. Kuster et al. [56] researched the influence of experimental drought and air warming on five years old saplings and found that bud break is advanced due to the carryover effect of the previous year’s drought.

In our study, there were eighteen examined phenological traits, describing duration and analog growth degree day sums for different periods important for bud flushing, as well as some of their ratios. Preliminary analysis showed no significant effect of trees within phenological groups, so variability of phenological traits was analyzed using two-way factorial analysis of variance, with Phenological group and Year as main effects. As expected, the effect of Phenological group was significant for all of them. The main effect of Year was not significant for any of the examined phenological traits, but the interaction Phenological group × Year was significant, suggesting considerable differences between two phenological groups in their reaction on conditions of examined years. In a rare multiannual study on the stability of phenological groups in pedunculate oak, Batos et al. [30] found preserved stability of leafing phenology in two populations near Belgrade despite differences in ecological conditions between three examined years, where 40%–50% of trees were classified in the same phenological group in all three years. However, they found significant effects of year for beginning, duration, and end of leafing in one population and only for duration in the other. Also, besides the fact that in general the two populations significantly differed in the beginning, duration, and end of leafing in all three years, in some cases that difference was not significant, suggesting differences in the reaction of populations on conditions in examined years. The best differentiation between interaction treatments was found in f01 trait, for which five homogenous groups were defined by Tukey’s test. Contributions to the total expected variance suggest that the strongest effect of phenological groups was found in f02a, f02b, and f02g, which supports their use in phenological variability studies. Also, strong interaction Phenological group × Year support, use of these characteristics in studies dealing with the phenological reaction of oak plants on differences between years. Two principal differences in reaction phenological groups on examined years were found, where years 2016 and 2017 had above average temperatures in March and April, respectively, while in 2015 the average monthly temperatures in both months were at the level of norm based on the period from 1961–1991 [46]. First, while in early group all phenological phases began earlier in warmer years 2016 and 2017 than in moderate 2015, in late group the only difference between 2015 and 2017 that was significant was for bud swelling (f01). Thus, earlier phenology in warmer years is clearer in the early bud-flushing group than in the late one. Second, examined periods between different phenological phases were significantly shorter in 2017 than in 2015 in early group, but significantly longer in the late group. According to the analysis of Batos et al. [30], in the vicinity of Belgrade the earlier beginning of the growing season in pedunculate oak is evident nowadays compared to the data from 40–50 years before. The same authors cited studies both in favor and against significant effect of years on phenology of oaks. Richardson et al. [19], in a detailed review paper, claimed that early phenological events (leafing, in particular) are closely related to climate change. They cite numerous works that promote early phenological events as indicators for tracking the effects of climate change even through centuries, as far as the tracking data spans, suggesting rapid rates of advance in temperate forests in recent decades (1.8–7.8 days/decade). However, for final conclusions on the effect of climate change on pedunculate oak early phenology, the research should be extended.

In Northern California, Herniman et al. [27] compared the phenology of four oak species—Blue oak (Quercus douglasii), California black oak (Q. kelloggii), Coast live oak (Q. agrifolia), and Oregon oak (Q. garryana) in 2014–2015 (a drought year) and 2015–2016 (“El Niño” year). They also compare phenophase onset date with growing degree days (GDD) and cumulative precipitation. The results show that in 2015 all four oak species had an earlier onset of breaking leaf buds, flowers and flower buds, and fruits compared to 2016. This is consistent with the research of Kuster et al. [56] on oak phenology, where the date of oak leaf bud burst has been shown to advance by 1–3 days per °C.

Trees require a specific amount of heat to develop from one point in their lifecycle to another. Studying genetic variation of the bud and leaf phenology in poplar clones during seven years, Pellis et al. [26] found that each year certain types of poplar clones opened buds earlier. They expressed bud burst as degree day sums (GDD) season. Each species requires a defined number of degree days to complete its development. In this research the patterns of genetic variation in leaf phenology have implications for short rotation intensive culture forestry and management of natural populations.

Considering their results, Knott et al. [28] found that 58-year-old populations of northern red oak (Quercus rubra) from USA Midwest, which originate from 32 locations across most of the native range in eastern North America areas with higher average springtime accumulated GDD, had significantly higher leaf-out sensitivity. Also, average springtime GDD is generally considered to be very relevant for predicting interannual variation in phenology.

In our study, GDD traits showed similar variability patterns as their analog duration traits. However, considerably higher effect of interaction Phenological group × Year on some phenological duration traits (f12b, f12g, f12a, f12a/f12g, and f12a/f12g) than on analog GDD traits, suggests that these duration traits would be preferred in studies concerning the difference in reaction of phenophorms on differences between years, particularly in studies concerning the effect of climate changes.

The majority of phenological traits showed poor broad sense heritability both in the early and late bud-flushing groups, suggesting poor differences within groups. Only some GDD traits, such as f02bt, f02gt, and f02at, achieved moderate values, which proposes them for further variability studies within phenophorm groups. Together with the clear difference between examined groups, these findings support the need for taxonomical differentiation between examined phenological groups, which should be further studied.

4.2. Relationship Between Examined Traits

In order to analyze relations between examined traits, the principal component analysis with rotation of selected principal components was performed. According to the highest loadings of original traits with selected and rotated principal components, six groups of traits were formed, with high correlation between traits within the same group and low correlation between traits that are not from the same group. In every group, there were traits from only one original class, except for acinx that was grouped in the first group with phenological traits, but its loading with the first principal component was moderate. Thus, phenological characteristics were in the first group, suggesting that they contribute the most to the total variance, but also, strong multicollinearity between them. Strong correlation between number of days and degree day sums needed for bud burst was also suggested by Pellis et al. [26]. This was expected because the study was performed in one area with similar conditions, but in populations on different areas this relationship could be different. Indeed, in the study of spatiotemporal pattern of Osmanthus fragrans phenology that examined five populations from South tropical to Warm temperate zone of China, Wang et al. [57] found that degree day sums in the certain period (from 1 January until 30 April) had strong negative correlation to period from 1 January until the date of bud burst and leafing phenophase. They also stressed the significance of cumulative precipitation and cumulative sunshine duration in the same period, as well as the strong effect of latitude. Leaf morphometric traits were distributed in the second (leaf size), fourth (leaf petiole based), and sixth (lbw/bl). All acorn traits were grouped in the third group, except for acinx in the first. The tree size characteristics were grouped in the fifth group, and there was no other trait included, suggesting that no other trait was in close relationship with tree size characteristics. It is also interesting that the only traits that had at least moderate loadings with the first principal component were acorn traits describing length (acl) and length/width ratio (acinx), suggesting that only these two traits were in moderate relationship with the group of examined phenological traits.

4.3. Relationship Between Trees

The analysis of the relationship between trees was based on 3D graph (Figure 4) formed by the first three principal components and dendrogram (Dendrogram 1) calculated by cluster analysis. Although it seems on 3D graph that all trees form one compact group, there is a clear distinction between early and late bud-flushing groups along the first principal component. All late trees are on the positive side of the PC1 axis, except for Tree 36, which is very close to zero, while all early trees are on the negative side. There is also considerable dispersion of trees by PC2 and PC3 axes within both early and late groups.

According to the hierarchical cluster analysis and scree test, five clusters were formed. All early trees were agglomerated in one cluster (Cluster number 2), together with late bud-flushing Tree 36, that was included in Cluster 2 probably due to leaf characteristics. All other clusters included trees from the late bud-flushing group, where the second cluster was the dominant cluster, while the other clusters consisted of just one tree.

Thus, the results of the study stress that distinct differentiation of examined trees according to the phenological groups suggest that examined phenological characteristics could have considerable importance in the description of the diversity of pedunculate oak.

However, both cluster analysis and principal component analysis indicate considerable contribution to the differences by leaf and acorn morphometric characteristics, too. Thus, all examined groups of traits could be considered informative in studies on the variability of pedunculate oak. The research should be continued on variability between populations and progenies, especially with respect to phenological and acorn morphometric traits. In light of climate change, the differential phenological responses of early and late oak genotypes could be key factors to consider when formulating adaptive forest management strategies.

5. Conclusions

According to the obtained results there was no clear distinction between early and late bud-burst groups of examined pedunculate oak trees in tree size and leaf morphometric characteristics, but it was in some acorn traits (acorn length and acorn index). The contribution of trees within Phenological groups to the total variance was the highest in leaf blade width (lbw) and leaf area (la). There was a clear effect Phenological group and interaction interaction Phenological group × Year but not of Years. It seems that the early bud flushing group reacts on warmer years with earlier phenology, which is not clear the late bud flushing group. Periods between different bud swelling and f12a, f12b, and f12g phenological phases were significantly shorter in 2017 than in 2015 in early group, but significantly longer in the late group. Principal component analysis revealed strict distinction between traits by their original nature, suggesting that all examined groups of traits could be considered as informative in studies on the variability of pedunculate oak. The tree size traits formed the separate, fifth group, suggesting no close relationship with any other examined characteristic. Both Cluster analysis and Principal component analysis suggest distinct classification by trees’ phenology, but also considerable differences by the second principal component which is closely related to leaf-size parameters. The research on variability between populations and progenies of pedunculate oak should be continued, especially in case of phenological and acorn morphometric traits. The differential phenological responses could be key factors in the design of adaptive forest management strategies, especially in the context of climate change.