European beech (Fagus sylvatica L.) seedlings (1–2 years old) were replanted from Denmark and Italy (Calabria). The established common garden experiment comprised three Italian provenances (distinct mountain ranges) with various altitudinal subgroups and two Danish provenances representing local populations. The seedlings were planted into HerkuPlast boxes—model Quick Pot 12 T/18 (plastic tray made by Kubern GmbH, Germany) in September 2016 and placed in an outdoor common garden experiment in the nursery of the University of Copenhagen. Plants from the same provenance subgroup were planted in each box (3×4 plants initially). Each block within the original common garden comprised at least nine Italian provenance subgroups (from three provenances), retaining an incomplete randomized block design (9 boxes in each of 20 incomplete blocks). The Italian provenances were replicated in 9 to 15 boxes. The Danish material was replicated within two blocks in 3 and 4 boxes of 12 seedlings each. Seedling mortality occurred after replanting, which violated the original design before the measurement date. The spatial design of the common garden experiment in September 2018 is visualized in Supplementary Figure S1. The imbalances due to this mortality were accounted for in the REML analysis (Butler et al., 2018), and the model details are specified in the description of the statistical analysis.

For our study, we selected Italian plant material from provenances Aspromonte and Sila Grande, each containing two within-provenance subgroups: Aspromonte (1,270 and 1,900 m.a.s.l.) and Sila Grande (1,300 and 1,900 m.a.s.l.) and Danish plant material originating from two different seed stands, F13—Stenderup Midtskov and F419—Holstenshuus (Table 1). The climate data in Table 1 are a selection of the bioclimatic variables from WorldClim version 2.1 climate data for 1970–2000 (Fick and Hijmans, 2017). All provenance origins and the experimental site are shown in Figure 1.

The main measurements (ChlF and GE) were done at the end of the second growing season of the trees in the common garden setting in September 2018 (10th–14th). There was sufficient precipitation at the end of August 2018 in Copenhagen, and during the measurement week, it was cloudy with occasional rainfall. The trees were not under any apparent stress during that time.

We measured ChlF kinetics using FluorPen FP110 (PSI, Brno, Czech Republic). We included all surviving (320) seedlings of the six provenance subgroups in the ChlF measurements (Table 2). ChlF induction (OJIP) parameters were measured after dark adaptation (detachable leaf clips placed on leaves 20 min before the measurement). After attaching the instrument to a leaf clip, a saturating light pulse (blue light, 470 nm) of 3,000 μmol·m⁻²·s⁻¹ intensity was applied to the leaf. Rapid fluorescence transient was captured within the first second of the measurement. We calculated the OJIP parameters according to Strasser et al. (2010) using the FluorPen software (Photon Systems Instruments, Brno, Czechia) and additional parameters from the I-P phase calculated using the raw data according to Oukarroum et al. (2009). All measured and calculated parameters are listed in Table 3.

We measured GE (representing net photosynthesis parameters) on randomly selected 84 individuals from all six provenance subgroups: A1270 (14), A1900 (17), SG1300 (19), SG1900 (19), DK-F13 (7), DK-F419 (8). We monitored the light curve using the gas exchange system LI-6400XT (LICOR, NE, United States), where a diode (6400-02B-LED) maintained uniform and stable light inside the chamber. The temperature was adjusted to the external conditions. We sampled leaves from the crown’s upper part for light curve measurement. The reference CO₂ concentration was set to 415 μmol·m⁻²·s⁻¹. We monitored the decreasing light as follows: 1,500-1,000-500-250-120-60-40-20-10-0 μmol·m⁻²·s⁻¹. The assimilation rate was displayed against the applied light, and the measured points were fitted according to the equation by Prado and De Moraes (1997) designed for woody species:

where P_N is the net photosynthesis rate [μmol (CO₂) m⁻²·s⁻¹]; P_gmax is the maximum gross photosynthesis rate [μmol (CO₂) m⁻²·s⁻¹] k is the adjusting factor; I is the photosynthetic photon flux density [μmol (photons) m⁻²·s⁻¹]; I_comp is the light compensation point [μmol (photons) m⁻²·s⁻¹]; and R_D is the dark respiration rate [μmol (CO₂) m⁻²·s⁻¹] derived from the curve according to Prado and De Moraes (1997). The constants were optimized by minimizing the residual sum of squares using the Solver function in MS Excel (Lobo et al., 2013). The GE parameters used in the analyses are listed in Table 4.

The first flushing assessment was done in May 2018 at the beginning of the second vegetation season. The second flushing assessment took place in April 2020 at a permanent site in a nearby forest where the seedlings were replanted in 2019. The flushing was scored only on the Italian provenances in 2018 and on all provenances in 2020. In total, 361 individuals were included in the first flushing assessment and 247 in the second. The flushing scores scale is based on a 7-point scale by Robson et al. (2013):

Leaf senescence was assessed at the beginning of November 2020 for all 244 individuals in the permanent experimental site. The senescence assessment was based on a 0–9 point scale assessing the majority of leaves: 0 represents leaves with no visible signs of senescence; 1–4 represents leaves with increasing levels of yellowing; 4–8 represents leaves with increasing amounts of necroses, and 9 represents dry and dead leaves.

The statistical analysis was conducted using the R software, version 4.2.1 (R Core Team, 2022).

A linear mixed model (LMM) was fitted to evaluate the photosynthetic parameters with the terms:

where Y corresponds to the data vector; μ is the overall mean effect; X is an incidence matrix of a fixed effect; a is the fixed vector of individual provenance; Z₁ and Z₂ are incidence matrices of random effects; S_i is the random effect of a block; S_i:S_j is the random effect of a box (S_j) nested in a block (S_i); e is the random vector of errors that were normally and independently distributed. The model was analyzed using the ASReml-R software (Butler et al., 2018). After fitting the model, outliers (0–10 per parameter) caused by error measurements were removed from the dataset. Custom R code for pairwise comparisons (least significant difference—LSD) was used for GE and ChlF parameters. We did not use any control for multiple testing.

Wald test (W-value) was implemented to evaluate the fixed effects of country of origin (a is the fixed vector of the country of origin in the model above for this analysis, replacing the fixed vector of individual provenance) when comparing their interaction on all measured traits between pooled Italian and Danish provenances.

For radar chart visualization, the following formula was used to perform a z-score normalization of the model-predicted mean values of photosynthetic parameters:

where X is the mean value of provenance in a given photosynthetic parameter predicted by the linear mixed model; μ is the mean of the data in a given photosynthetic parameter; σ is the standard deviation of the data in a given photosynthetic parameter.

The flushing scores were analyzed using a two-sided Wilcoxon signed-rank test 4.2.1 (R Core Team, 2022).

The Italian provenances had significantly higher average values compared to the Danish ones in ChlF parameters F_m/F₀ (W = 1.55e-06), F_v/F_m (W = 3.93e-07), ψ₀ (W = 1.53e-05), φ_E0 (W = 7.27e-08), δ_R0 (W = 0.001), φ_R0 (W = 2.76–07), PI_ABS (W = 5.16e-08), PI_TOTAL (W = 6.58e-07), ΔV_IP (W = 3.75e-06), and significantly lower average values than the Danish ones in parameters φ_D0 (W = 3.76e-07), ABS/RC (W = 4.24e-07), TR₀/RC (W = 1.47e-05), DI₀/RC (W = 8.81e-09). There was no clear observable difference in parameters ET₀/RC (W = 0.721), φ_Pav (W = 0.259). The model-predicted means of ChlF parameters for each provenance subgroup are outlined in Table 5. The Italian provenances manifested higher average values than the Danish ones in GE parameters LSP₉₅ (W = 7.90e-08) and P_N (W = 2.64e-05). Under the same irradiation, Italian provenances reached the maximum assimilation rate faster than the Danish ones, and their net assimilation was higher. There was no clear trend or difference in the light compensation point (LCP) (W = 0.764) or dark respiration (R_D) (W = 0.309). All provenances needed the same photon flux density to assimilate more CO₂ than they respire. No difference was observed in the quantum yield of photosynthesis (ϕ_I0) (W = 0.701), which signifies that the molar ratio between carbon assimilated to photons absorbed in the process does not differ among provenances. The model-predicted means of GE parameters for each provenance subgroup are shown in Table 6. The performance of individual provenance subgroups for each chosen parameter from both types of measurement is visualized in Figure 2.

All provenance subgroup pairs were compared in photosynthetic parameters of ChlF and GE (Figure 3). Both ChlF and GE parameters distinguished Italian and Danish provenances. ChlF parameters DI0/RC, φ_R0, PI_TOTAL, and ΔV_IP differed among Italian provenance subgroups based on altitude (Figure 3–pink). The GE parameters LSP₉₅ and P_N showed differentiation between Italian provenances from Aspromonte and Sila Grande, i.e., different geographical locations (Figure 3–gray and red). In contrast, no altitudinal trend was visible from the GE data, while ChlF measurements could not separate all geographically distant provenances. A different pattern occurred when comparing the two local Danish provenances. The ChlF measurements yielded seven usable parameters (F_m/F₀, F_v/F_m, φ_E0, φ_D0, ABS/RC, TR₀/RC, DI₀/RC) that capture differences between the two Danish provenances (Figure 3–light green). The two Danish provenances could not be separated based on the GE parameters. The only parameter that distinguished the Italian within-provenance variation in altitude of origin and also the Danish provenances was DI₀/RC (Figure 3–pink and light green).

In May 2018, provenances from the two mountain ranges, Sila Grande (SG) and Aspromonte (A), varied significantly (p = 1.92E-15) in bud flushing scoring. The more southern Aspromonte showed more developed buds at the time of assessment. An altitudinal difference could be seen within the Sila Grande provenance, where SG1900 showed significantly less (p = 2.87E-14) developed buds than SG1300. No such difference could be observed within the Aspromonte provenance. The differences from May 2018 between all individual pairs are outlined in Table 7. The distribution of flushing scores illustrating the magnitude of difference between the provenance subgroups in May 2018 is shown in Figure 4A. The flushing assessment from April 2020 on the surviving trees revealed the same persisting trend (visualized in Figure 4B), although the flushing was in an earlier stage during this assessment. The differences between the provenances and their subgroups in April 2020 (Table 7) remained significant for the same pairs except for SG1300 and A1270. The overall difference between Aspromonte and Sila Grande remained substantial (Figure 4–red vs. green). In addition, in 2020, the Danish provenances were also scored. Both Danish provenances were in a later flushing stage than the Italian ones, except for SG1300, which showed a late flushing pattern more similar to the Danish beech. The significance of the differences is shown in Table 7.

At the beginning of November 2020, all the provenances showed signs of senescence. The majority of leaves were yellowing on almost all of the assessed trees. No significant differences were captured between the provenance subgroups. The senescence scores for individual provenance subgroups are visualized in Figure 5.

We highlighted the significant correlations (α = 0.05) parameters obtained from ChlF and GE methods in a heat map (Supplementary Figure S2). There were mostly weak correlations between ChlF and GE parameters, with the highest positive being 0.41 between PI_TOTAL and P_N and the highest negative being −0.40 between TR₀/RC and P_N.

Other studies comparing the photosynthetic performance of European beech that originates from different altitudes include provenances where the trend of elevation is confounded in a shared cardinal direction (e.g., Robson et al., 2012, 2013; Kučerová et al., 2018; Pšidová et al., 2018). Therefore, the latitude and longitude of origin could also drive the differences and cannot be effectively separated. In our study, the geographical distances between provenance subgroups occupying a single mountain range were minimal, meaning that physiological differences within these provenances reflect mainly altitude.

The most common use of ChlF kinetics is to diagnose early stress responses to biotic or abiotic factors before any decline in growth can be observed (Fiorani and Schurr, 2013; Faseela et al., 2020). DI₀/RC and F_v/Fm are often used as drought or temperature stress indicators (Mathur et al., 2013; Yuan et al., 2013; Fghire et al., 2015). However, the most commonly used parameter, F_v/F_m, is quite robust and sometimes cannot be as sensitive as other ChlF parameters (Bussotti et al., 2020). As the trees in our experiment were under no apparent stress at the time of measurements, these parameters instead reflected provenance variation and overall vitality after replanting. When comparing the subgroups of Italian provenances in ChlF, four parameters could make a distinction based on the altitude (DI₀/RC, φ_R0, PI_TOTAL, ΔV_IP). Higher φ_R0, PI_TOTAL, and ΔV_IP values were estimated for the high-altitude provenance subgroups in Aspromonte and Sila Grande, and higher DI₀/RC values for the low-altitude provenance subgroups in Aspromonte and Sila Grande (Table 5). ΔV_IP is the parameter we found especially interesting in which the relative loss of the I-P phase is related to a loss of PSI reaction centers (Oukarroum et al., 2009), as indices connected to the I-P phase were not much previously explored in European beech. Furthermore, F_v/F_m and ΔV_IP are suggested to be the two most suitable OJIP parameters to capture the variability of plant photosynthetic efficiency and their responses to environmental pressures (Bussotti et al., 2020). The ChlF parameter we found most indicative was DI₀/RC, as it was the only one to differentiate between the two Danish provenances as well as between the Italian provenance subgroups on an altitudinal basis. However, in our study, ChlF parameters could not distinguish provenances based on their geographical origin.

Findings from two studies conducted in a single common garden experiment in Slovakia by Kučerová et al. (2018) and Pšidová et al. (2018) reported mechanistic differences in ChlF kinetics among beech provenances. Kučerová et al. (2018) observed that Austrian and Polish provenances differed from others in the F_v/F_m ratio, although this trend was inconsistent along the altitudinal gradient. Pšidová et al. (2018), focused on provenance performance during heatwaves, finding that ABS/RC and ψ₀ were more limited in lower altitude provenances. Higher altitude provenances showed better PSII photochemistry performance, indicated by PI_ABS, under both dry and wet conditions, which aligns with our findings. Kurjak et al. (2019) found significant associations between ChlF parameters (mainly F_v/F_m and PI_ABS) and the climate of origin. Provenances from colder, wetter regions exhibited higher PSII performance compared to those from warmer, drier climates. This trend was also linked to the geographical distance from the Slovenian refugium and the common garden site, with closer proximity resulting in better PSII performance. PSII performance decreased towards the north-western margins, supporting our results regarding latitude-related variations in ChlF.

On the other hand, two GE parameters (LSP₉₅ and P_N) could reliably distinguish geographical origin (the different Italian mountain ranges, Aspromonte and Sila Grande). We observed no trend in altitude in the GE data. In contrast, Kučerová et al. (2018) reported a significant linear altitudinal trend for the CO₂ assimilation rate between Central European provenances of European beech. Although the provenances used in Kučerová et al. (2018) study originated from different altitudes, they were also geographically distant; therefore, the latitude and longitude of origin could also drive the differences. The GE measurements revealed higher LSP₉₅ and P_N values for more southern Aspromonte than Sila Grande. The higher LSP₉₅ values from Aspromonte suggest that Sila Grande beech may be more shade tolerant. A study by Meng et al. (2014) showed that shade-adapted plants have lower light saturation levels. However, in beech, the release of sunlight after a windstorm tends to diminish the differences in LSP between fully released (formerly shade-grown) and sun-grown seedlings (Reynolds et al., 2003). In our study, the different LSP₉₅ values of Aspromonte and Sila Grande could most likely be attributed to the local adaptation. For example, Čater and Levanič (2019) report a reduced shade tolerance and higher efficiency of beech towards the south, which supports our results and may correspond to a boost in regeneration.

In 2018 and 2020, the southern provenance from Aspromonte exhibited significantly higher flushing scores compared to Sila Grande’s provenance. Additionally, in 2020, the Danish provenances showed lower flushing scores than the Italian ones. Previous studies have found that flushing is significantly correlated with latitude, longitude, and altitude, influenced by factors such as chilling, temperature sums, and day length (Nielsen and Jørgensen, 2003). The duration of phenological phases varies annually and mainly depends on temperature (Bednářová and Merklová, 2007; Slovíková and Bednářová, 2014). Similarly, research from Slovakia (Schieber et al., 2013), southern Germany (Dittmar and Elling, 2006), and Slovenia (Čufar et al., 2012) showed phenological variation in beech along the altitudinal gradient. In southern Germany, beech leaf unfolding was delayed by 2–3 days for every 100 m of altitudinal increase (Dittmar and Elling, 2006). While the Sila Grande provenance showed variation in flushing along the elevation gradient in 2018 and 2020, this was not observed within Aspromonte. However, our consistent results from 2018 and 2020 suggest that the adaptation of the provenances is linked to flushing variation, as the provenance flushing scores remained stable in Danish climatic conditions over two seasons. Vitasse et al. (2009) noted that beech leaf unfolding is less temperature-sensitive than other European tree species, indicating that other factors are more influential.

Flushing scores did not correlate significantly with photosynthetic parameters. The two data collections in 2018 (flushing vs. GE and ChlF) were conducted at the beginning and end of the vegetation season, respectively, so the expected correlation pattern may have been disrupted due to the various factors influencing the trees during the vegetation season. The weather fluctuations influenced the seedlings’ performance in the common garden, and some mortality occurred during the vegetation season. However, both GE and flushing similarly cluster the Italian provenances and their subgroups (Figures 2, 4). In 2020, the Danish provenances were also scored. Both Danish provenances were in a later flushing stage than the Italian ones. The only exception was SG1300, which showed a late flushing pattern similar to the Danish beech. It still had a slightly earlier flushing pattern than the Danish provenances, although not significantly.

There were no visible differences in the senescence of leaves between the provenances and provenance subgroups at the beginning of November 2020 (Figure 5). We conclude that as there were no observable differences in November senescence, the photosynthetic measurements conducted at the beginning of September 2018 were not influenced by the direct onset of leaf senescence. It is widely known that senescence is connected to chlorophyll content. As leaves in deciduous trees lose some chlorophyll content during senescence, the photosynthetic capacity decreases on a leaf area basis (Adams et al., 1990). Leaves that fall later in the season have a better translocation of nutrients, and the translocation is species-dependent (Niinemets and Tamm, 2005). A study in another beech species (Fagus grandifolia Ehrh.) documented chlorophyll degradation that was not sudden but continuous, as the slower formation of the abscission layer during senescence allowed beech to translocate more nutrients out of leaves (Primka and Smith, 2019). However, ChlF does not seem to be affected to such a degree by large decreases in chlorophyll content as the yellowing leaves of deciduous trees are still efficiently used in photosynthesis because the photosynthetic capacity on a chlorophyll basis and photosynthetic energy conversion efficiency do not show a decline (Adams et al., 1990). Therefore, we believe that even at the beginning of the senescence process, the ChlF method could still yield relevant results.

This study has potential limitations that need to be stated. Firstly, we only used two provenances from Calabria representing beech refugia, which we separated into two distinct altitudinal subgroups and two local Danish provenances. Secondly, the sample size of GE measurements was smaller than that of ChlF measurements due to the extensive time consumption of the GE method. The main measurements (GE and ChlF) were conducted during one week at the beginning of September, while flushing had to be scored in the Spring (2018 and 2020). Our study did not involve any stress treatment or biomass measurements, which would also be beneficial in understanding the well-being of replanted trees in new conditions. Flushing was not scored for the Danish plant material in 2018. Senescence was only assessed in 2020; however, this secondary evaluation provided more insight into the potential influence of different onsets of senescence on the measurements taken in September 2018. Our results and conclusions are based on methods applied at various times, as we aimed to observe the differences in spring and autumn phenology. This combined approach helped to verify our main conclusion connected to the variance of photosynthetic performance due to phenology.